1 DATA

1.1 Pop sive over time dataset

data <- read.table(file=here::here("data", "BE_Census_Data.csv"), sep=",", header=TRUE)
  
head(data)
##    Block       Old_Label                  Label Population Genetic_Diversity
## 1 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High1      High1              High
## 2 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High2      High2              High
## 3 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High3      High3              High
## 4 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High4      High4              High
## 5 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High5      High5              High
## 6 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High6      High6              High
##   Generation_Parents Generation_Eggs Date_Census Date_Start_Eggs Date_End_Eggs
## 1                  0               1     2/17/21         2/17/21       2/18/21
## 2                  0               1     2/17/21         2/17/21       2/18/21
## 3                  0               1     2/17/21         2/17/21       2/18/21
## 4                  0               1     2/17/21         2/17/21       2/18/21
## 5                  0               1     2/17/21         2/17/21       2/18/21
## 6                  0               1     2/17/21         2/17/21       2/18/21
##   Image_Box1 Image_Box2 MC_Box1 MC_Box2 MC_Total_Adults HC_Box1 HC_Box2
## 1   DSC_0155       <NA>      NA      NA           101.9      NA      NA
## 2   DSC_0156       <NA>      NA      NA           104.7      NA      NA
## 3   DSC_0157       <NA>      NA      NA           105.8      NA      NA
## 4   DSC_0158       <NA>      NA      NA           101.9      NA      NA
## 5   DSC_0159       <NA>      NA      NA           102.2      NA      NA
## 6   DSC_0160       <NA>      NA      NA           101.7      NA      NA
##   HC_Total_Adults Nb_parents Growth_rate First_Throw Extinction_finalstatus
## 1             100         NA          NA          no                     no
## 2             100         NA          NA          no                     no
## 3             100         NA          NA          no                     no
## 4             100         NA          NA          no                     no
## 5             100         NA          NA          no                     no
## 6             100         NA          NA          no                     no
##   Less_than_5
## 1          no
## 2          no
## 3          no
## 4          no
## 5          no
## 6          no
tail(data)
##      Block         Old_Label                  Label Population
## 571 Block5 B5 - Backup Fam H 8/25 BE B5 | G6 Med 19      Med19
## 572 Block5 B5 - Backup Fam I 8/25 BE B5 | G6 Med 20      Med20
## 573 Block5   B5-Backup Fam A 8/25 BE B5 | G6 Med 21      Med21
## 574 Block5   B5-Backup Fam B 8/25 BE B5 | G6 Med 22      Med22
## 575 Block5   B5-Backup Fam C 8/25 BE B5 | G6 Med 23      Med23
## 576 Block5   B5-Backup Fam D 8/25 BE B5 | G6 Med 24      Med24
##     Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 571               Med                  5               6     8/25/21
## 572               Med                  5               6     8/25/21
## 573               Med                  5               6     8/25/21
## 574               Med                  5               6     8/25/21
## 575               Med                  5               6     8/25/21
## 576               Med                  5               6     8/25/21
##     Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2  MC_Box1  MC_Box2
## 571         8/25/21       8/26/21   DSC_0273   DSC_0274  2.00000  2.00000
## 572         8/25/21       8/26/21       <NA>       <NA>       NA       NA
## 573         8/25/21       8/26/21   DSC_0277   DSC_0278 29.76978 32.92378
## 574         8/25/21       8/26/21       <NA>       <NA>       NA       NA
## 575         8/25/21       8/26/21   DSC_0267   DSC_0268 25.31021 46.10092
## 576         8/25/21       8/26/21   DSC_0283   DSC_0284 71.09684 55.37443
##     MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 571             4.0       2       2               4          8    0.500000
## 572             0.0      NA      NA              NA          0          NA
## 573            62.7      23      23              46         35    1.314286
## 574             0.0      NA      NA              NA         NA          NA
## 575            71.4      20      36              56         40    1.400000
## 576           126.5      76      57             133         41    3.243902
##     First_Throw Extinction_finalstatus Less_than_5
## 571          no                     no         yes
## 572     extinct                    yes         yes
## 573          no                     no          no
## 574     extinct                    yes         yes
## 575          no                     no          no
## 576          no                     no          no
dim(data)
## [1] 576  23
summary(data)
##     Block            Old_Label            Label            Population       
##  Length:576         Length:576         Length:576         Length:576        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##  Genetic_Diversity  Generation_Parents Generation_Eggs Date_Census       
##  Length:576         Min.   :0.0        Min.   :1.0     Length:576        
##  Class :character   1st Qu.:1.0        1st Qu.:2.0     Class :character  
##  Mode  :character   Median :2.5        Median :3.5     Mode  :character  
##                     Mean   :2.5        Mean   :3.5                       
##                     3rd Qu.:4.0        3rd Qu.:5.0                       
##                     Max.   :5.0        Max.   :6.0                       
##                                                                          
##  Date_Start_Eggs    Date_End_Eggs       Image_Box1         Image_Box2       
##  Length:576         Length:576         Length:576         Length:576        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##     MC_Box1          MC_Box2       MC_Total_Adults     HC_Box1      
##  Min.   :  0.00   Min.   :  0.00   Min.   :  0.00   Min.   :  0.00  
##  1st Qu.: 26.52   1st Qu.: 25.99   1st Qu.: 51.05   1st Qu.: 22.00  
##  Median : 58.87   Median : 54.29   Median :101.50   Median : 56.50  
##  Mean   : 75.88   Mean   : 70.57   Mean   :128.37   Mean   : 72.06  
##  3rd Qu.:106.93   3rd Qu.: 99.40   3rd Qu.:170.25   3rd Qu.:101.00  
##  Max.   :327.40   Max.   :299.00   Max.   :514.40   Max.   :288.00  
##  NA's   :142      NA's   :141      NA's   :5        NA's   :144     
##     HC_Box2       HC_Total_Adults   Nb_parents      Growth_rate    
##  Min.   :  0.00   Min.   :  0.0   Min.   :  0.00   Min.   :0.0000  
##  1st Qu.: 25.00   1st Qu.: 62.0   1st Qu.: 67.25   1st Qu.:0.4183  
##  Median : 52.00   Median :100.0   Median :100.00   Median :0.9386  
##  Mean   : 68.08   Mean   :132.1   Mean   :138.30   Mean   :1.4482  
##  3rd Qu.: 96.75   3rd Qu.:169.0   3rd Qu.:188.00   3rd Qu.:2.3587  
##  Max.   :290.00   Max.   :499.0   Max.   :499.00   Max.   :6.8571  
##  NA's   :142      NA's   :44      NA's   :122      NA's   :142     
##  First_Throw        Extinction_finalstatus Less_than_5       
##  Length:576         Length:576             Length:576        
##  Class :character   Class :character       Class :character  
##  Mode  :character   Mode  :character       Mode  :character  
##                                                              
##                                                              
##                                                              
## 
#Remove populations killed due to the pathogen
data_status  <- data.frame(data$Population,data$Extinction_finalstatus)
data_kill<-data[data$Extinction_finalstatus=="kill",]
dim(data_kill)
## [1] 72 23
data<-data[data$Extinction_finalstatus!="kill",]

#Add pop growth rate
data$Nb_adults <- data$HC_Total_Adults
data$Lambda <- data$Nb_adults / data$Nb_parents
data$obs<-as.factor(1:nrow(data))
#summary(data)
#Remove 

#Check variables
str(data)
## 'data.frame':    504 obs. of  26 variables:
##  $ Block                 : chr  "Block3" "Block3" "Block3" "Block3" ...
##  $ Old_Label             : chr  "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " ...
##  $ Label                 : chr  "BE B3 2/17 | G1  High1" "BE B3 2/17 | G1  High2" "BE B3 2/17 | G1  High3" "BE B3 2/17 | G1  High4" ...
##  $ Population            : chr  "High1" "High2" "High3" "High4" ...
##  $ Genetic_Diversity     : chr  "High" "High" "High" "High" ...
##  $ Generation_Parents    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Generation_Eggs       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Date_Census           : chr  "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
##  $ Date_Start_Eggs       : chr  "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
##  $ Date_End_Eggs         : chr  "2/18/21" "2/18/21" "2/18/21" "2/18/21" ...
##  $ Image_Box1            : chr  "DSC_0155" "DSC_0156" "DSC_0157" "DSC_0158" ...
##  $ Image_Box2            : chr  NA NA NA NA ...
##  $ MC_Box1               : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ MC_Box2               : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ MC_Total_Adults       : num  102 105 106 102 102 ...
##  $ HC_Box1               : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ HC_Box2               : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ HC_Total_Adults       : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ Nb_parents            : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Growth_rate           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ First_Throw           : chr  "no" "no" "no" "no" ...
##  $ Extinction_finalstatus: chr  "no" "no" "no" "no" ...
##  $ Less_than_5           : chr  "no" "no" "no" "no" ...
##  $ Nb_adults             : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ Lambda                : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ obs                   : Factor w/ 504 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
data$Genetic_Diversity <- as.factor(data$Genetic_Diversity)
data$Block <- as.factor(data$Block)
data$Population <- as.factor(data$Population)
data$Extinction_finalstatus <- as.factor(data$Extinction_finalstatus)


#Order levels 
data$Genetic_Diversity <- plyr::revalue(data$Genetic_Diversity, c("Med"="Medium"))
data$Genetic_Diversity <- factor(data$Genetic_Diversity, levels = c("High", "Medium", "Low"))
levels(data$Genetic_Diversity)
## [1] "High"   "Medium" "Low"
data<-droplevels(data) #remove previous levels 
str(data)
## 'data.frame':    504 obs. of  26 variables:
##  $ Block                 : Factor w/ 3 levels "Block3","Block4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Old_Label             : chr  "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " ...
##  $ Label                 : chr  "BE B3 2/17 | G1  High1" "BE B3 2/17 | G1  High2" "BE B3 2/17 | G1  High3" "BE B3 2/17 | G1  High4" ...
##  $ Population            : Factor w/ 84 levels "High1","High13",..: 1 9 19 24 25 26 27 28 39 49 ...
##  $ Genetic_Diversity     : Factor w/ 3 levels "High","Medium",..: 1 1 1 1 1 1 1 3 3 3 ...
##  $ Generation_Parents    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Generation_Eggs       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Date_Census           : chr  "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
##  $ Date_Start_Eggs       : chr  "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
##  $ Date_End_Eggs         : chr  "2/18/21" "2/18/21" "2/18/21" "2/18/21" ...
##  $ Image_Box1            : chr  "DSC_0155" "DSC_0156" "DSC_0157" "DSC_0158" ...
##  $ Image_Box2            : chr  NA NA NA NA ...
##  $ MC_Box1               : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ MC_Box2               : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ MC_Total_Adults       : num  102 105 106 102 102 ...
##  $ HC_Box1               : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ HC_Box2               : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ HC_Total_Adults       : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ Nb_parents            : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Growth_rate           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ First_Throw           : chr  "no" "no" "no" "no" ...
##  $ Extinction_finalstatus: Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Less_than_5           : chr  "no" "no" "no" "no" ...
##  $ Nb_adults             : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ Lambda                : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ obs                   : Factor w/ 504 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
head(data)
##    Block       Old_Label                  Label Population Genetic_Diversity
## 1 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High1      High1              High
## 2 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High2      High2              High
## 3 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High3      High3              High
## 4 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High4      High4              High
## 5 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High5      High5              High
## 6 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High6      High6              High
##   Generation_Parents Generation_Eggs Date_Census Date_Start_Eggs Date_End_Eggs
## 1                  0               1     2/17/21         2/17/21       2/18/21
## 2                  0               1     2/17/21         2/17/21       2/18/21
## 3                  0               1     2/17/21         2/17/21       2/18/21
## 4                  0               1     2/17/21         2/17/21       2/18/21
## 5                  0               1     2/17/21         2/17/21       2/18/21
## 6                  0               1     2/17/21         2/17/21       2/18/21
##   Image_Box1 Image_Box2 MC_Box1 MC_Box2 MC_Total_Adults HC_Box1 HC_Box2
## 1   DSC_0155       <NA>      NA      NA           101.9      NA      NA
## 2   DSC_0156       <NA>      NA      NA           104.7      NA      NA
## 3   DSC_0157       <NA>      NA      NA           105.8      NA      NA
## 4   DSC_0158       <NA>      NA      NA           101.9      NA      NA
## 5   DSC_0159       <NA>      NA      NA           102.2      NA      NA
## 6   DSC_0160       <NA>      NA      NA           101.7      NA      NA
##   HC_Total_Adults Nb_parents Growth_rate First_Throw Extinction_finalstatus
## 1             100         NA          NA          no                     no
## 2             100         NA          NA          no                     no
## 3             100         NA          NA          no                     no
## 4             100         NA          NA          no                     no
## 5             100         NA          NA          no                     no
## 6             100         NA          NA          no                     no
##   Less_than_5 Nb_adults Lambda obs
## 1          no       100     NA   1
## 2          no       100     NA   2
## 3          no       100     NA   3
## 4          no       100     NA   4
## 5          no       100     NA   5
## 6          no       100     NA   6
tail(data)
##      Block         Old_Label                  Label Population
## 571 Block5 B5 - Backup Fam H 8/25 BE B5 | G6 Med 19      Med19
## 572 Block5 B5 - Backup Fam I 8/25 BE B5 | G6 Med 20      Med20
## 573 Block5   B5-Backup Fam A 8/25 BE B5 | G6 Med 21      Med21
## 574 Block5   B5-Backup Fam B 8/25 BE B5 | G6 Med 22      Med22
## 575 Block5   B5-Backup Fam C 8/25 BE B5 | G6 Med 23      Med23
## 576 Block5   B5-Backup Fam D 8/25 BE B5 | G6 Med 24      Med24
##     Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 571            Medium                  5               6     8/25/21
## 572            Medium                  5               6     8/25/21
## 573            Medium                  5               6     8/25/21
## 574            Medium                  5               6     8/25/21
## 575            Medium                  5               6     8/25/21
## 576            Medium                  5               6     8/25/21
##     Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2  MC_Box1  MC_Box2
## 571         8/25/21       8/26/21   DSC_0273   DSC_0274  2.00000  2.00000
## 572         8/25/21       8/26/21       <NA>       <NA>       NA       NA
## 573         8/25/21       8/26/21   DSC_0277   DSC_0278 29.76978 32.92378
## 574         8/25/21       8/26/21       <NA>       <NA>       NA       NA
## 575         8/25/21       8/26/21   DSC_0267   DSC_0268 25.31021 46.10092
## 576         8/25/21       8/26/21   DSC_0283   DSC_0284 71.09684 55.37443
##     MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 571             4.0       2       2               4          8    0.500000
## 572             0.0      NA      NA              NA          0          NA
## 573            62.7      23      23              46         35    1.314286
## 574             0.0      NA      NA              NA         NA          NA
## 575            71.4      20      36              56         40    1.400000
## 576           126.5      76      57             133         41    3.243902
##     First_Throw Extinction_finalstatus Less_than_5 Nb_adults   Lambda obs
## 571          no                     no         yes         4 0.500000 499
## 572     extinct                    yes         yes        NA       NA 500
## 573          no                     no          no        46 1.314286 501
## 574     extinct                    yes         yes        NA       NA 502
## 575          no                     no          no        56 1.400000 503
## 576          no                     no          no       133 3.243902 504

1.2 Remaining heterozygosity

#Upload growth rate end of experiment
data_he <- read.table(file=here::here("data", "He_ReMaining_BeatricePop.csv"), sep=",", header=TRUE)
data_he$Population <- as.factor(data_he$Population)
data_he <- data_he[,c(1,3)]


### Additional: He remaining 
# New vector for the lost during exp
data_he$He_remain_exp <- NA

# He lost during exp
for(i in levels(data$Population)){
temp_data_pop <- subset(data, Population==i)
ifelse(sum(temp_data_pop$Generation_Parents)!=15,print("ERROR"),print(i))
temp_he_remaining_gen <- 1 - (1/(2*temp_data_pop$Nb_parents))

#To consider extinct add this row: 
is.na(temp_he_remaining_gen) <- sapply (temp_he_remaining_gen, is.infinite)

temp_he_remain_exp <- prod(temp_he_remaining_gen, na.rm = TRUE)
data_he$He_remain_exp[data_he$Population==i] <- temp_he_remain_exp
rm(temp_he_remain_exp,temp_he_remaining_gen,temp_data_pop)
}
## [1] "High1"
## [1] "High13"
## [1] "High14"
## [1] "High15"
## [1] "High16"
## [1] "High17"
## [1] "High18"
## [1] "High19"
## [1] "High2"
## [1] "High20"
## [1] "High21"
## [1] "High22"
## [1] "High23"
## [1] "High24"
## [1] "High25"
## [1] "High26"
## [1] "High27"
## [1] "High28"
## [1] "High3"
## [1] "High32"
## [1] "High33"
## [1] "High34"
## [1] "High35"
## [1] "High4"
## [1] "High5"
## [1] "High6"
## [1] "High7"
## [1] "Low1"
## [1] "Low10"
## [1] "Low11"
## [1] "Low12"
## [1] "Low13"
## [1] "Low14"
## [1] "Low15"
## [1] "Low16"
## [1] "Low17"
## [1] "Low18"
## [1] "Low19"
## [1] "Low2"
## [1] "Low20"
## [1] "Low21"
## [1] "Low22"
## [1] "Low24"
## [1] "Low25"
## [1] "Low26"
## [1] "Low27"
## [1] "Low28"
## [1] "Low29"
## [1] "Low3"
## [1] "Low30"
## [1] "Low31"
## [1] "Low32"
## [1] "Low33"
## [1] "Low34"
## [1] "Low35"
## [1] "Low36"
## [1] "Low4"
## [1] "Low5"
## [1] "Low7"
## [1] "Low8"
## [1] "Low9"
## [1] "Med1"
## [1] "Med10"
## [1] "Med11"
## [1] "Med12"
## [1] "Med13"
## [1] "Med14"
## [1] "Med15"
## [1] "Med17"
## [1] "Med18"
## [1] "Med19"
## [1] "Med2"
## [1] "Med20"
## [1] "Med21"
## [1] "Med22"
## [1] "Med23"
## [1] "Med24"
## [1] "Med3"
## [1] "Med4"
## [1] "Med5"
## [1] "Med6"
## [1] "Med7"
## [1] "Med8"
## [1] "Med9"
# Remaining He at the end of the experiment
data_he$He_remain_end <- data_he$He_remain * data_he$He_remain_exp

# Remove kill population 
data_he <- data_he[!is.na(data_he$He_remain_exp),]
data_he <- droplevels(data_he)

1.3 Phenotyping dataset

#Upload growth rate end of experiment
data_phenotyping <- read.table(file=here::here("data", "Adaptation_Data.csv"), sep=",", header=TRUE)


#Factors variables 
data_phenotyping$ID_Rep <- as.factor(data_phenotyping$ID_Rep)
data_phenotyping$Week <- as.factor(data_phenotyping$Week)
data_phenotyping$Block <- as.factor(data_phenotyping$Block)
data_phenotyping$Population <- as.factor(data_phenotyping$Population)
#Revalue
data_phenotyping$Genetic_Diversity <- as.factor(data_phenotyping$Divsersity)
data_phenotyping$Genetic_Diversity <-   plyr::revalue(data_phenotyping$Genetic_Diversity, c("Med"="Medium"))
data_phenotyping$Genetic_Diversity <- factor(data_phenotyping$Genetic_Diversity, levels = c("High", "Medium", "Low"))



#Add sum total 
data_phenotyping$Nb_ttx <- rowSums(data_phenotyping[,11:13], na.rm=TRUE)
#Proportion
data_phenotyping$p_adults <- data_phenotyping$Adults  / data_phenotyping$Nb_ttx
data_phenotyping$p_pupae <- data_phenotyping$Pupae / data_phenotyping$Nb_ttx
data_phenotyping$p_larvae <- data_phenotyping$Larvae / data_phenotyping$Nb_ttx

data_phenotyping_all <- data_phenotyping
data_phenotyping <- data_phenotyping[data_phenotyping$Week=="5-week",]

1.4 Merge dataset He

#Add heterozygosity
data <- merge(x = data, y = data_he, by = "Population", all.x=TRUE)

#Add heterozygosity

#To add population extinct
data_phenotyping_withextinct <- merge(x = data_phenotyping, y = data_he, by = "Population", all.y=TRUE)
data_phenotyping <- merge(x = data_phenotyping, y = data_he, by = "Population", all.x=TRUE)

1.5 Subset dataset generation

#Data beginning
data_G2 <- data[data$Generation_Eggs==2,]
#Data end 
data_G6 <- data[data$Generation_Eggs==6,]

1.6 Merge dataset phenotyping

## Merge dataset phenotyping
#merge dataset
temp_Start <- data_G2[,c("Block","Population",
                         "Genetic_Diversity","Lambda", "He_remain", "He_remain_end" )]
temp_Start$ID_Rep <- 1
temp_Start$Phenotyping <- "Initial"

temp_end <- data_phenotyping[,c("Block","Population",
                                "Genetic_Diversity","Lambda","ID_Rep", "He_remain", "He_remain_end" )]
temp_end$Phenotyping <- "Final"

data_evolution_Lambda <- rbind(temp_Start,temp_end)
data_evolution_Lambda$Phenotyping <- factor(data_evolution_Lambda$Phenotyping, 
                                            levels = c("Initial", "Final"))

2 PLOT

2.1 Population size

PLOT_Pop_size<-ggplot2::ggplot(data, aes(x = factor(Generation_Eggs), 
                                y = Nb_adults, 
                                group = Population,
                                colour = Genetic_Diversity)) + 
  geom_point(position = position_dodge(0.4), size = 0.8, alpha = 0.7) +
  geom_line(position = position_dodge(0.4), size = 0.2, alpha = 0.6) +
  ylab("Population size") +
  labs(color="Genetic diversity") +
  xlab("Generation") +
  theme_LO
PLOT_Pop_size

# cowplot::save_plot(file = here::here("figures","OldPlot_PopulationSize.pdf"),PLOT_Pop_size, base_height = 10/cm(1),
#           base_width = 15/cm(1), dpi = 610)

SUM_Popsize <- Rmisc::summarySE(data, 
                        measurevar = c("Nb_adults"),
                       groupvars = c("Generation_Eggs",
                                     "Genetic_Diversity"), 
                       na.rm=TRUE)
SUM_Popsize
##    Generation_Eggs Genetic_Diversity  N Nb_adults        sd        se       ci
## 1                1              High 27 100.00000   0.00000  0.000000  0.00000
## 2                1            Medium 23 100.00000   0.00000  0.000000  0.00000
## 3                1               Low 34 100.00000   0.00000  0.000000  0.00000
## 4                2              High 27 344.29630  73.77294 14.197610 29.18360
## 5                2            Medium 23 282.04348 100.75735 21.009360 43.57075
## 6                2               Low 34 282.61765  78.53006 13.467794 27.40043
## 7                3              High 27 151.85185  80.96899 15.582489 32.03027
## 8                3            Medium 23 118.73913  88.54491 18.462891 38.28969
## 9                3               Low 34 104.61765  85.13341 14.600260 29.70445
## 10               4              High 26 114.00000  80.20224 15.728954 32.39439
## 11               4            Medium 23  62.65217  64.50483 13.450188 27.89398
## 12               4               Low 34  46.38235  39.63241  6.796903 13.82840
## 13               5              High 27 103.62963  51.56486  9.923661 20.39838
## 14               5            Medium 20  61.50000  50.59904 11.314290 23.68108
## 15               5               Low 31  51.51613  46.13377  8.285870 16.92200
## 16               6              High 27 147.66667  42.60282  8.198916 16.85311
## 17               6            Medium 19  69.73684  46.02396 10.558620 22.18284
## 18               6               Low 28  77.03571  46.14428  8.720450 17.89289
# SUM_Popsize <- Rmisc::summarySE(data[data$Extinction_finalstatus=="no",],
#                         measurevar = c("Nb_adults"),
#                        groupvars = c("Generation_Eggs",
#                                      "Genetic_Diversity"),
#                        na.rm=TRUE)
# SUM_Popsize



PLOT_Pop_size_average <- PLOT_Pop_size + 
  geom_point(data = SUM_Popsize, aes(x = factor(Generation_Eggs), 
                                y = Nb_adults, 
                                group = Genetic_Diversity,
                                colour = Genetic_Diversity), 
             size = 5, position = position_dodge(0.4)) +
  geom_errorbar(data = SUM_Popsize, aes(x = factor(Generation_Eggs), 
                                ymin = Nb_adults-ci, ymax = Nb_adults+ci, 
                                group = Genetic_Diversity,
                                colour = Genetic_Diversity), 
             size = 1, width=.2,  position = position_dodge(0.4)) +
  geom_line(data = SUM_Popsize, aes(x = factor(Generation_Eggs), 
                                y = Nb_adults, 
                                group = Genetic_Diversity,
                                colour = Genetic_Diversity), 
           linetype = "dashed", size = 1.5, position = position_dodge(0.4)) 
PLOT_Pop_size_average

# cowplot::save_plot(file = here::here("figures","PopulationSize_average.pdf"),   PLOT_Pop_size_average, base_height = 10/cm(1),           base_width = 15/cm(1), dpi = 610)

2.2 Extinction

## Extinction yes no
PLOT_Extinction<-ggplot2::ggplot(data, aes(x = factor(Generation_Eggs), 
                                y = Nb_adults, 
                                group = Population,
                                colour = Extinction_finalstatus)) + 
  facet_wrap(~Genetic_Diversity, ncol=3) + 
  #geom_point(position = position_dodge(0.4), size = 0.4, alpha = 0.7) +
  geom_line(position = position_dodge(0.1), size = 0.4, alpha = 0.8) +
  ylab("Population size") +
  scale_color_manual(values = c("black", "red")) +
  xlab("Generation") +
  theme(legend.position = "none") +
  theme_LO
PLOT_Extinction

# #
# cowplot::save_plot(here::here("figures","Extinction.pdf"),   PLOT_Extinction, base_height = 6/cm(1),
#           base_width = 20/cm(1), dpi = 610)
# 
# #

2.3 Growth rate G2

tapply(data_G2$Lambda,data_G2$Genetic_Diversity,mean)
##     High   Medium      Low 
## 3.442963 2.820435 2.826176
PLOT_Growth<-ggplot2::ggplot(data_G2, aes(x = Genetic_Diversity, y = Lambda, 
                         colour = Genetic_Diversity)) + 
  geom_boxplot(aes(group = Genetic_Diversity),outlier.colour = NA) +
  geom_jitter(width = 0.25, size = 1, alpha = 0.8) +
  ylab("Growth rate") +
  theme(legend.position = "none") +
  xlab("Genetic diversity") +
  theme_LO
PLOT_Growth

# 
# cowplot::save_plot(here::here("figures","OldPlot_Growth_Start.pdf"),   PLOT_Growth, base_height = 10/cm(1),
#           base_width = 8/cm(1), dpi = 610)

2.4 Growth rate evolution

SUM_MEAN_start_end <- Rmisc::summarySE(data_evolution_Lambda, measurevar = c("Lambda"),
                       groupvars = c("Genetic_Diversity","Phenotyping"), 
                       na.rm=TRUE)
SUM_POP_start_end <- Rmisc::summarySE(data_evolution_Lambda, measurevar = c("Lambda"),
                       groupvars = c("Genetic_Diversity","Phenotyping",
                                     "Population"), 
                       na.rm=TRUE)





PLOT_Growth_startend<-ggplot(SUM_POP_start_end, aes(x = Phenotyping, 
                                 y = Lambda, 
                                 group = Population,
                                 colour = Genetic_Diversity)) + 
  geom_point(position = position_dodge(0.4), size = 0.8, alpha = 0.7) +
  geom_line(position = position_dodge(0.4), size = 0.2, alpha = 0.6) +
  ylab("Growth rate") +
  xlab("Phenotyping") +
  theme_LO
PLOT_Growth_startend

PLOT_Growth_startend_average <- PLOT_Growth_startend + 
  geom_point(data = SUM_MEAN_start_end, aes(x = Phenotyping, 
                                y = Lambda, 
                                group = Genetic_Diversity,
                                colour = Genetic_Diversity), 
             size = 5, position = position_dodge(0.4)) +
  geom_errorbar(data = SUM_MEAN_start_end, aes(x = Phenotyping, 
                                ymin = Lambda-ci, ymax = Lambda+ci, 
                                group = Genetic_Diversity,
                                colour = Genetic_Diversity), 
             size = 1, width=.2,  position = position_dodge(0.4)) +
  geom_line(data = SUM_MEAN_start_end, aes(x = Phenotyping, 
                                y = Lambda, 
                                group = Genetic_Diversity,
                                colour = Genetic_Diversity), 
           linetype = "dashed", size = 1.5, position = position_dodge(0.4)) 
PLOT_Growth_startend_average

# 
# cowplot::save_plot(here::here("figures","OldPlot_Lambda_Evolution.pdf"),
#           PLOT_Growth_startend_average, base_height = 10/cm(1),
#           base_width = 15/cm(1), dpi = 610)

2.5 Life stage proportion

SUM_Prop_adults<- Rmisc::summarySE(data_phenotyping_all, 
                        measurevar = c("p_adults"),
                       groupvars = c("Genetic_Diversity","Week"), 
                       na.rm = TRUE)
SUM_Prop_adults
##   Genetic_Diversity   Week   N  p_adults         sd          se         ci
## 1              High 4-week  50 0.6508098 0.17776618 0.025139934 0.05052059
## 2              High 5-week 105 0.9632040 0.05259208 0.005132461 0.01017786
## 3            Medium 4-week  24 0.7121270 0.15722597 0.032093616 0.06639070
## 4            Medium 5-week  49 0.9563136 0.06160975 0.008801393 0.01769639
## 5               Low 4-week  30 0.6568606 0.19029917 0.034743716 0.07105888
## 6               Low 5-week  59 0.9582021 0.07845154 0.010213520 0.02044458
SUM_Prop_pupae<- Rmisc::summarySE(data_phenotyping_all, 
                        measurevar = c("p_pupae"),
                       groupvars = c("Genetic_Diversity","Week"), 
                       na.rm = TRUE)
SUM_Prop_pupae
##   Genetic_Diversity   Week   N    p_pupae         sd          se          ci
## 1              High 4-week  50 0.28825555 0.15629878 0.022103985 0.044419621
## 2              High 5-week 105 0.01458769 0.02407445 0.002349426 0.004658999
## 3            Medium 4-week  24 0.23019148 0.14506173 0.029610601 0.061254195
## 4            Medium 5-week  49 0.02283167 0.03485130 0.004978757 0.010010462
## 5               Low 4-week  30 0.29380227 0.17562943 0.032065401 0.065581108
## 6               Low 5-week  59 0.01675308 0.02459640 0.003202178 0.006409856
SUM_Prop_larvae<- Rmisc::summarySE(data_phenotyping_all, 
                        measurevar = c("p_larvae"),
                       groupvars = c("Genetic_Diversity","Week"), 
                       na.rm = TRUE)
SUM_Prop_larvae
##   Genetic_Diversity   Week   N   p_larvae         sd          se          ci
## 1              High 4-week  50 0.06093466 0.04274850 0.006045551 0.012148989
## 2              High 5-week 105 0.02220828 0.04324043 0.004219833 0.008368088
## 3            Medium 4-week  24 0.05768150 0.05999020 0.012245449 0.025331640
## 4            Medium 5-week  49 0.02085474 0.04462923 0.006375605 0.012819012
## 5               Low 4-week  30 0.04933713 0.08198954 0.014969174 0.030615398
## 6               Low 5-week  59 0.02504479 0.06691319 0.008711355 0.017437671
SUM_Prop <- data.frame(SUM_Prop_adults[,1:4],
                       p_pupae=SUM_Prop_pupae$p_pupae,
                       p_larvae=SUM_Prop_larvae$p_larvae)


SUM_Prop<-tidyr::gather(SUM_Prop,"Stage", "Proportion",4:6,-"Genetic_Diversity")
SUM_Prop
##    Genetic_Diversity   Week   N    Stage Proportion
## 1               High 4-week  50 p_adults 0.65080979
## 2               High 5-week 105 p_adults 0.96320403
## 3             Medium 4-week  24 p_adults 0.71212702
## 4             Medium 5-week  49 p_adults 0.95631359
## 5                Low 4-week  30 p_adults 0.65686059
## 6                Low 5-week  59 p_adults 0.95820213
## 7               High 4-week  50  p_pupae 0.28825555
## 8               High 5-week 105  p_pupae 0.01458769
## 9             Medium 4-week  24  p_pupae 0.23019148
## 10            Medium 5-week  49  p_pupae 0.02283167
## 11               Low 4-week  30  p_pupae 0.29380227
## 12               Low 5-week  59  p_pupae 0.01675308
## 13              High 4-week  50 p_larvae 0.06093466
## 14              High 5-week 105 p_larvae 0.02220828
## 15            Medium 4-week  24 p_larvae 0.05768150
## 16            Medium 5-week  49 p_larvae 0.02085474
## 17               Low 4-week  30 p_larvae 0.04933713
## 18               Low 5-week  59 p_larvae 0.02504479
SUM_Prop$Stage <- factor(SUM_Prop$Stage, levels = c("p_adults", "p_pupae", "p_larvae"))


PLOT_prop <- ggplot2::ggplot(data = SUM_Prop, aes(x = Genetic_Diversity, 
                                                  y = Proportion, fill = Stage)) +
  facet_wrap(~ Week, ncol=2) +
  geom_bar(stat="identity") + 
  xlab("Level of genetic diversity") +
  labs(color="Stage of individuals") +
  ylab("Proportion of each stage") +
  scale_fill_manual(values= c("#4EABB9", "#E95C2B", "#CBA938"), 
                    breaks = c("p_adults", "p_pupae", "p_larvae"),
                    labels = c( "Adult", "Pupae","Larvae")) + 
  theme_LO
PLOT_prop

# 
# cowplot::save_plot(here::here("figures","Life_stage_proportion.pdf"),   PLOT_prop,
#           base_height = 8/cm(1), base_width = 16/cm(1), dpi = 610)

2.6 Remaining heterozygosity

# All but with: He beginning for both and with pseudoreplication for final
ggplot2::ggplot(data_evolution_Lambda, aes(x = He_remain, y = Lambda, 
                         colour = Genetic_Diversity)) + 
  facet_wrap(~ Phenotyping, ncol = 1, scales = "free")+
  geom_point(size = 2, alpha = 0.8) +
  ylab("Growth rate") +
  theme(legend.position = "none") +
  xlab("Remaining heterozygosity") +
  theme_LO

PLOT_He_Initial <- ggplot2::ggplot(data_G2, aes(x = He_remain, y = Lambda, 
                         colour = Genetic_Diversity)) + 
  geom_point(size = 2, alpha = 0.8) +
  ylab("Growth rate") +
  theme(legend.position = "none") +
  xlab("Remaining heterozygosity") +
  theme_LO
PLOT_He_Initial

PLOT_He_End <- ggplot2::ggplot(data_evolution_Lambda[data_evolution_Lambda$Phenotyping=="Final",],
                               aes(x = He_remain_end, y = Lambda, 
                         colour = Genetic_Diversity)) + 
  geom_point(size = 2, alpha = 0.8) +
  ylab("Growth rate") +
  theme(legend.position = "none") +
  xlab("Remaining heterozygosity") +
  theme_LO
PLOT_He_End

# Remove pseudoreplication: average of each population
SUM_POP_He_Mean<- Rmisc::summarySE(data_evolution_Lambda[data_evolution_Lambda$Phenotyping=="Final",],
                            measurevar = c("Lambda"),
                       groupvars = c("Genetic_Diversity","He_remain_end",
                                     "Population"), 
                       na.rm=TRUE)

PLOT_He_Final <- ggplot2::ggplot(SUM_POP_He_Mean, aes(x = He_remain_end, y = Lambda, 
                         colour = Genetic_Diversity)) + 
  geom_point(size = 2, alpha = 0.8) +
  ylab("Growth rate") +
  theme(legend.position = "none") +
  xlab("Remaining heterozygosity") +
  theme_LO
PLOT_He_Final

# ALL WITHOUT PSEUDO
SUM_POP_He_ALL<- Rmisc::summarySE(data_evolution_Lambda,
                            measurevar = c("Lambda"),
                       groupvars = c("Genetic_Diversity","Phenotyping", "He_remain", "He_remain_end",
                                     "Population"), 
                       na.rm=TRUE)

SUM_POP_He_ALL$HE <- ifelse(SUM_POP_He_ALL$Phenotyping=="Initial",SUM_POP_He_ALL$He_remain,SUM_POP_He_ALL$He_remain_end)

SUM_POP_He_ALL <- SUM_POP_He_ALL[!is.na(SUM_POP_He_ALL$He_remain_end),]
SUM_POP_He_ALL <- SUM_POP_He_ALL[SUM_POP_He_ALL$HE!="-Inf",]


facet_names <- c('Initial'="Before rescue experiment",'Final'="After rescue experiment")

PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL, aes(x = HE, y = Lambda, 
                         colour = Genetic_Diversity, 
                         shape = Phenotyping)) +
  #facet_wrap(~ Phenotyping, ncol=1) +
  geom_point(size = 3, alpha = 0.8) +
  scale_shape_manual(values = c(1, 19)) + 
  ylab("Growth rate") +
  xlab("Remaining heterozygosity") +
  theme_LO
PLOT_He_ALL

PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL, aes(x = HE, y = Lambda, 
                         colour = Genetic_Diversity, 
                         shape = Phenotyping)) +
  facet_wrap(~ Phenotyping, ncol=1) +
  geom_point(size = 3, alpha = 0.8) +
  scale_shape_manual(values = c(1, 19)) + 
  ylab("Growth rate") +
  xlab("Remaining heterozygosity") +
  theme_LO
PLOT_He_ALL

PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL, aes(x = HE, y = Lambda,
                         colour = Genetic_Diversity)) +
  facet_wrap(~ Phenotyping, ncol=1) +
  geom_point(size = 3, alpha = 0.8) +
  ylab("Growth rate") +
  xlab("Remaining heterozygosity") +
  theme_LO
PLOT_He_ALL

# cowplot::save_plot(here::here("figures","Fitness_He_initial_final.pdf"),   PLOT_He_ALL,
#           base_height = 20/cm(1), base_width = 16/cm(1), dpi = 610)



# PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",], aes(x = HE, y = Lambda,
#                          colour = Genetic_Diversity,
#                          shape = Phenotyping)) +
#   geom_point(size = 3, alpha = 0.8) +
#   scale_shape_manual(values = c(19)) +
#   ylab("Growth rate") +
#   xlab("Remaining heterozygosity") +
#   theme_LO
# PLOT_He_ALL



########################
###Plot final phenotyping with extinct 
########################
data_phenotyping_withextinct
##     Population   Week  Block ID_Rep Divsersity Population.1 Box Initial.N
## 1        High1 5-week Block3     52       High            1   1        30
## 2        High1 5-week Block3     53       High            1   2        30
## 3       High13 5-week Block4    130       High           13   2        30
## 4       High13 5-week Block4    131       High           13   3        30
## 5       High13 5-week Block4    129       High           13   1        30
## 6       High14 5-week Block4    132       High           14   1        30
## 7       High14 5-week Block4    133       High           14   2        30
## 8       High14 5-week Block4    134       High           14   3        30
## 9       High14 5-week Block4    135       High           14   4        30
## 10      High15 5-week Block4    136       High           15   1        30
## 11      High15 5-week Block4    137       High           15   2        30
## 12      High16 5-week Block4    138       High           16   1        30
## 13      High16 5-week Block4    139       High           16   2        30
## 14      High16 5-week Block4    140       High           16   3        30
## 15      High16 5-week Block4    141       High           16   4        30
## 16      High16 5-week Block4    142       High           16   5        30
## 17      High17 5-week Block4    143       High           17   1        30
## 18      High18 5-week Block4    145       High           18   2        30
## 19      High18 5-week Block4    146       High           18   3        30
## 20      High18 5-week Block4    144       High           18   1        30
## 21      High19 5-week Block4    147       High           19   1        30
## 22      High19 5-week Block4    149       High           19   3        30
## 23      High19 5-week Block4    148       High           19   2        30
## 24       High2 5-week Block3     54       High            2   1        30
## 25       High2 5-week Block3     55       High            2   2        30
## 26      High20 5-week Block4    150       High           20   1        30
## 27      High20 5-week Block4    151       High           20   2        30
## 28      High20 5-week Block4    152       High           20   3        30
## 29      High20 5-week Block4    153       High           20   4        30
## 30      High20 5-week Block4    154       High           20   5        30
## 31      High21 5-week Block4    155       High           21   1        30
## 32      High21 5-week Block4    157       High           21   3        30
## 33      High21 5-week Block4    156       High           21   2        30
## 34      High22 5-week Block4    158       High           22   1        30
## 35      High22 5-week Block4    159       High           22   2        30
## 36      High23 5-week Block4    160       High           23   1        30
## 37      High23 5-week Block4    162       High           23   3        30
## 38      High23 5-week Block4    163       High           23   4        30
## 39      High23 5-week Block4    161       High           23   2        30
## 40      High24 5-week Block4    164       High           24   1        30
## 41      High24 5-week Block4    165       High           24   2        30
## 42      High25 5-week Block5    185       High           25   3        30
## 43      High25 5-week Block5    184       High           25   2        30
## 44      High25 5-week Block5    183       High           25   1        30
## 45      High25 5-week Block5    186       High           25   4        30
## 46      High25 5-week Block5    187       High           25   5        30
## 47      High26 5-week Block5    188       High           26   1        30
## 48      High27 5-week Block5    190       High           27   2        30
## 49      High27 5-week Block5    189       High           27   1        30
## 50      High27 5-week Block5    191       High           27   3        30
## 51      High27 5-week Block5    192       High           27   4        30
## 52      High28 5-week Block5    198       High           28   6        30
## 53      High28 5-week Block5    193       High           28   1        30
## 54      High28 5-week Block5    196       High           28   4        30
## 55      High28 5-week Block5    195       High           28   3        30
## 56      High28 5-week Block5    194       High           28   2        30
## 57      High28 5-week Block5    197       High           28   5        30
## 58       High3 5-week Block3     56       High            3   1        30
## 59       High3 5-week Block3     57       High            3   2        30
## 60       High3 5-week Block3     58       High            3   3        30
## 61      High32 5-week Block5    203       High           32   1        30
## 62      High32 5-week Block5    204       High           32   2        30
## 63      High33 5-week Block5    205       High           33   1        30
## 64      High33 5-week Block5    206       High           33   2        30
## 65      High34 5-week Block5    210       High           34   4        30
## 66      High34 5-week Block5    211       High           34   5        30
## 67      High34 5-week Block5    209       High           34   3        30
## 68      High34 5-week Block5    207       High           34   1        30
## 69      High34 5-week Block5    208       High           34   2        30
## 70      High35 5-week Block5    215       High           35   4        30
## 71      High35 5-week Block5    212       High           35   1        30
## 72      High35 5-week Block5    213       High           35   2        30
## 73      High35 5-week Block5    214       High           35   3        30
## 74      High35 5-week Block5    217       High           35   6        30
## 75      High35 5-week Block5    216       High           35   5        30
## 76       High4 5-week Block3     60       High            4   2        30
## 77       High4 5-week Block3     61       High            4   3        30
## 78       High4 5-week Block3     62       High            4   4        30
## 79       High4 5-week Block3     59       High            4   1        30
## 80       High4 5-week Block3     63       High            4   5        30
## 81       High5 5-week Block3     64       High            5   1        30
## 82       High5 5-week Block3     65       High            5   2        30
## 83       High5 5-week Block3     66       High            5   3        30
## 84       High5 5-week Block3     67       High            5   4        30
## 85       High5 5-week Block3     68       High            5   5        30
## 86       High5 5-week Block3     69       High            5   6        30
## 87       High5 5-week Block3     70       High            5   7        30
## 88       High6 5-week Block3     72       High            6   2        30
## 89       High6 5-week Block3     73       High            6   3        30
## 90       High6 5-week Block3     71       High            6   1        30
## 91       High7 5-week Block3     74       High            7   1        30
## 92       High7 5-week Block3     75       High            7   2        30
## 93       High7 5-week Block3     77       High            7   4        30
## 94       High7 5-week Block3     78       High            7   5        30
## 95       High7 5-week Block3     79       High            7   6        30
## 96       High7 5-week Block3     76       High            7   3        30
## 97        Low1 5-week Block3      5        Low            1   5        30
## 98        Low1 5-week Block3      1        Low            1   1        30
## 99        Low1 5-week Block3      2        Low            1   2        30
## 100       Low1 5-week Block3      3        Low            1   3        30
## 101       Low1 5-week Block3      4        Low            1   4        30
## 102      Low10 5-week Block3     29        Low           10   1        30
## 103      Low11 5-week Block3     30        Low           11   1        30
## 104      Low12 5-week Block3     32        Low           12   2        30
## 105      Low12 5-week Block3     31        Low           12   1        30
## 106      Low13 5-week Block4     85        Low           13   1        30
## 107      Low14 5-week Block4     87        Low           14   2        30
## 108      Low14 5-week Block4     88        Low           14   3        30
## 109      Low14 5-week Block4     86        Low           14   1        30
## 110      Low15 5-week Block4     90        Low           15   2        30
## 111      Low15 5-week Block4     89        Low           15   1        30
## 112      Low16 5-week Block4     91        Low           16   1        30
## 113      Low17 5-week Block4     93        Low           17   2        30
## 114      Low17 5-week Block4     92        Low           17   1        30
## 115      Low18 5-week Block4     97        Low           18   4        30
## 116      Low18 5-week Block4     94        Low           18   1        30
## 117      Low18 5-week Block4     95        Low           18   2        30
## 118      Low18 5-week Block4     96        Low           18   3        30
## 119      Low19 5-week Block4     98        Low           19   1        30
## 120       Low2 5-week Block3      7        Low            2   2        30
## 121       Low2 5-week Block3      8        Low            2   3        30
## 122       Low2 5-week Block3      9        Low            2   4        30
## 123       Low2 5-week Block3      6        Low            2   1        30
## 124      Low20 5-week Block4     99        Low           20   1        30
## 125      Low21 5-week Block4    100        Low           21   1        30
## 126      Low22 5-week Block4    101        Low           22   1        30
## 127      Low24 5-week Block4    103        Low           24   1        30
## 128      Low25 5-week Block4    105        Low           25   2        30
## 129      Low25 5-week Block4    104        Low           25   1        30
## 130      Low26 5-week Block5    166        Low           26   1        30
## 131      Low29 5-week Block5    167        Low           29   1        30
## 132       Low3 5-week Block3     11        Low            3   2        30
## 133       Low3 5-week Block3     12        Low            3   3        30
## 134       Low3 5-week Block3     10        Low            3   1        30
## 135      Low33 5-week Block5    168        Low           33   1        30
## 136      Low34 5-week Block5    171        Low           34   3        30
## 137      Low34 5-week Block5    172        Low           34   4        30
## 138      Low34 5-week Block5    170        Low           34   2        30
## 139      Low34 5-week Block5    169        Low           34   1        30
## 140      Low35 5-week Block5    173        Low           35   1        30
## 141       Low4 5-week Block3     13        Low            4   1        30
## 142       Low4 5-week Block3     14        Low            4   2        30
## 143       Low4 5-week Block3     15        Low            4   3        30
## 144       Low4 5-week Block3     16        Low            4   4        30
## 145       Low4 5-week Block3     17        Low            4   5        30
## 146       Low5 5-week Block3     19        Low            5   2        30
## 147       Low5 5-week Block3     20        Low            5   3        30
## 148       Low5 5-week Block3     18        Low            5   1        30
## 149       Low8 5-week Block3     21        Low            8   1        30
## 150       Low8 5-week Block3     22        Low            8   2        30
## 151       Low8 5-week Block3     24        Low            8   4        30
## 152       Low8 5-week Block3     23        Low            8   3        30
## 153       Low9 5-week Block3     25        Low            9   1        30
## 154       Low9 5-week Block3     26        Low            9   2        30
## 155       Low9 5-week Block3     27        Low            9   3        30
## 156       Low9 5-week Block3     28        Low            9   4        30
## 157       Med1 5-week Block3     33        Med            1   1        30
## 158      Med10 5-week Block4    110        Med           10   1        30
## 159      Med11 5-week Block4    112        Med           11   2        30
## 160      Med11 5-week Block4    113        Med           11   3        30
## 161      Med11 5-week Block4    111        Med           11   1        30
## 162      Med11 5-week Block4    114        Med           11   4        30
## 163      Med12 5-week Block4    115        Med           12   1        30
## 164      Med13 5-week Block4    118        Med           13   3        30
## 165      Med13 5-week Block4    116        Med           13   1        30
## 166      Med13 5-week Block4    117        Med           13   2        30
## 167      Med15 5-week Block4    119        Med           15   1        30
## 168      Med15 5-week Block4    122        Med           15   4        30
## 169      Med15 5-week Block4    123        Med           15   5        30
## 170      Med15 5-week Block4    120        Med           15   2        30
## 171      Med15 5-week Block4    121        Med           15   3        30
## 172      Med18 5-week Block5    174        Med           18   1        30
## 173      Med18 5-week Block5    176        Med           18   3        30
## 174      Med18 5-week Block5    177        Med           18   4        30
## 175      Med18 5-week Block5    175        Med           18   2        30
## 176      Med19 5-week Block5    178        Med           19   1        30
## 177       Med2 5-week Block3     34        Med            2   1        30
## 178       Med2 5-week Block3     35        Med            2   2        30
## 179       Med2 5-week Block3     36        Med            2   3        30
## 180      Med21 5-week Block5    179        Med           21   1        30
## 181      Med23 5-week Block5    180        Med           23   1        30
## 182      Med24 5-week Block5    181        Med           24   1        30
## 183      Med24 5-week Block5    182        Med           24   2        30
## 184       Med3 5-week Block3     37        Med            3   1        30
## 185       Med3 5-week Block3     38        Med            3   2        30
## 186       Med3 5-week Block3     39        Med            3   3        30
## 187       Med4 5-week Block3     40        Med            4   1        30
## 188       Med5 5-week Block3     41        Med            5   1        30
## 189       Med6 5-week Block3     44        Med            6   3        30
## 190       Med6 5-week Block3     42        Med            6   1        30
## 191       Med6 5-week Block3     43        Med            6   2        30
## 192       Med7 5-week Block3     46        Med            7   2        30
## 193       Med7 5-week Block3     45        Med            7   1        30
## 194       Med8 5-week Block3     47        Med            8   1        30
## 195       Med8 5-week Block3     48        Med            8   2        30
## 196       Med8 5-week Block3     49        Med            8   3        30
## 197       Med8 5-week Block3     51        Med            8   5        30
## 198       Med8 5-week Block3     50        Med            8   4        30
## 199       Med9 5-week Block4    106        Med            9   1        30
## 200       Med9 5-week Block4    109        Med            9   4        30
## 201       Med9 5-week Block4    107        Med            9   2        30
## 202       Med9 5-week Block4    108        Med            9   3        30
## 203      Low27   <NA>   <NA>   <NA>       <NA>           NA  NA        NA
## 204      Low28   <NA>   <NA>   <NA>       <NA>           NA  NA        NA
## 205      Low30   <NA>   <NA>   <NA>       <NA>           NA  NA        NA
## 206      Low31   <NA>   <NA>   <NA>       <NA>           NA  NA        NA
## 207      Low32   <NA>   <NA>   <NA>       <NA>           NA  NA        NA
## 208      Low36   <NA>   <NA>   <NA>       <NA>           NA  NA        NA
## 209       Low7   <NA>   <NA>   <NA>       <NA>           NA  NA        NA
## 210      Med14   <NA>   <NA>   <NA>       <NA>           NA  NA        NA
## 211      Med17   <NA>   <NA>   <NA>       <NA>           NA  NA        NA
## 212      Med20   <NA>   <NA>   <NA>       <NA>           NA  NA        NA
## 213      Med22   <NA>   <NA>   <NA>       <NA>           NA  NA        NA
##       Start     End Larvae Pupae Adults     Lambda Genetic_Diversity Nb_ttx
## 1    7/8/21 8/12/21      3     9     97 3.23333333              High    109
## 2    7/8/21 8/12/21      1     1      8 0.26666667              High     10
## 3   7/15/21 8/19/21      1     2    104 3.46666667              High    107
## 4   7/15/21 8/19/21      2     1     84 2.80000000              High     87
## 5   7/15/21 8/19/21      2     0     84 2.80000000              High     86
## 6   7/15/21 8/19/21      1     0     83 2.76666667              High     84
## 7   7/15/21 8/19/21      2     0     48 1.60000000              High     50
## 8   7/15/21 8/19/21      0     0     62 2.06666667              High     62
## 9   7/15/21 8/19/21      1     0     43 1.43333333              High     44
## 10  7/15/21 8/19/21      4     0     84 2.80000000              High     88
## 11  7/15/21 8/19/21      8     0     57 1.90000000              High     65
## 12  7/15/21 8/19/21      1     1    110 3.66666667              High    112
## 13  7/15/21 8/19/21      1     0     45 1.50000000              High     46
## 14  7/15/21 8/19/21      0     0     59 1.96666667              High     59
## 15  7/15/21 8/19/21      3     0     54 1.80000000              High     57
## 16  7/15/21 8/19/21      1     2     56 1.86666667              High     59
## 17  7/15/21 8/19/21      0     3    165 5.50000000              High    168
## 18  7/15/21 8/19/21      3     2     50 1.66666667              High     55
## 19  7/15/21 8/19/21      0     0     88 2.93333333              High     88
## 20  7/15/21 8/19/21      0     5     75 2.50000000              High     80
## 21  7/15/21 8/19/21      3     6     66 2.20000000              High     75
## 22  7/15/21 8/19/21      0     1     93 3.10000000              High     94
## 23  7/15/21 8/19/21      1     1    112 3.73333333              High    114
## 24   7/8/21 8/12/21      4     0    165 5.50000000              High    169
## 25   7/8/21 8/12/21      0     1     73 2.43333333              High     74
## 26  7/15/21 8/19/21      5     0     77 2.56666667              High     82
## 27  7/15/21 8/19/21      1     1     84 2.80000000              High     86
## 28  7/15/21 8/19/21      0     0    101 3.36666667              High    101
## 29  7/15/21 8/19/21      2     1     48 1.60000000              High     51
## 30  7/15/21 8/19/21      0     0     80 2.66666667              High     80
## 31  7/15/21 8/19/21      0     0     92 3.06666667              High     92
## 32  7/15/21 8/19/21      1     4     86 2.86666667              High     91
## 33  7/15/21 8/19/21      2     2    121 4.03333333              High    125
## 34  7/15/21 8/19/21      2     4    135 4.50000000              High    141
## 35  7/15/21 8/19/21      0     0     83 2.76666667              High     83
## 36  7/15/21 8/19/21      1     0     61 2.03333333              High     62
## 37  7/15/21 8/19/21      0     0     86 2.86666667              High     86
## 38  7/15/21 8/19/21      0     0     47 1.56666667              High     47
## 39  7/15/21 8/19/21      0     0     66 2.20000000              High     66
## 40  7/15/21 8/19/21      2     5     88 2.93333333              High     95
## 41  7/15/21 8/19/21     12     3    105 3.50000000              High    120
## 42  7/22/21 8/26/21      2     1     98 3.26666667              High    101
## 43  7/22/21 8/26/21      1     5     71 2.36666667              High     77
## 44  7/22/21 8/26/21      2     1     59 1.96666667              High     62
## 45  7/22/21 8/26/21      4     2     86 2.86666667              High     92
## 46  7/22/21 8/26/21      1     0     91 3.03333333              High     92
## 47  7/22/21 8/26/21      0     0     79 2.63333333              High     79
## 48  7/22/21 8/26/21      4     3     71 2.36666667              High     78
## 49  7/22/21 8/26/21     37     4     73 2.43333333              High    114
## 50  7/22/21 8/26/21      1     2     73 2.43333333              High     76
## 51  7/22/21 8/26/21      2     7     83 2.76666667              High     92
## 52  7/22/21 8/26/21      0     0     23 0.76666667              High     23
## 53  7/22/21 8/26/21      1     1    103 3.43333333              High    105
## 54  7/22/21 8/26/21      1     3     65 2.16666667              High     69
## 55  7/22/21 8/26/21      1     0     75 2.50000000              High     76
## 56  7/22/21 8/26/21      1     0    129 4.30000000              High    130
## 57  7/22/21 8/26/21      0     0     56 1.86666667              High     56
## 58   7/8/21 8/12/21      1     0    124 4.13333333              High    125
## 59   7/8/21 8/12/21     16     0     72 2.40000000              High     88
## 60   7/8/21 8/12/21      0     0     35 1.16666667              High     35
## 61  7/22/21 8/26/21      3     2    136 4.53333333              High    141
## 62  7/22/21 8/26/21      2     0    143 4.76666667              High    145
## 63  7/22/21 8/26/21      1     2     61 2.03333333              High     64
## 64  7/22/21 8/26/21      2     0     36 1.20000000              High     38
## 65  7/22/21 8/26/21      0     0    116 3.86666667              High    116
## 66  7/22/21 8/26/21      0     0     36 1.20000000              High     36
## 67  7/22/21 8/26/21      0     1     70 2.33333333              High     71
## 68  7/22/21 8/26/21      1     1    101 3.36666667              High    103
## 69  7/22/21 8/26/21      0     0     82 2.73333333              High     82
## 70  7/22/21 8/26/21      1     0    101 3.36666667              High    102
## 71  7/22/21 8/26/21      0     0    107 3.56666667              High    107
## 72  7/22/21 8/26/21      0     0     49 1.63333333              High     49
## 73  7/22/21 8/26/21      1     0     75 2.50000000              High     76
## 74  7/22/21 8/26/21      0     0     36 1.20000000              High     36
## 75  7/22/21 8/26/21      0     0     48 1.60000000              High     48
## 76   7/8/21 8/12/21      0     0    186 6.20000000              High    186
## 77   7/8/21 8/12/21      0     0    107 3.56666667              High    107
## 78   7/8/21 8/12/21      0     0     73 2.43333333              High     73
## 79   7/8/21 8/12/21      0     0    145 4.83333333              High    145
## 80   7/8/21 8/12/21      0     0     35 1.16666667              High     35
## 81   7/8/21 8/12/21      0     1     44 1.46666667              High     45
## 82   7/8/21 8/12/21      0     0     93 3.10000000              High     93
## 83   7/8/21 8/12/21      0     0     84 2.80000000              High     84
## 84   7/8/21 8/12/21      1     0     81 2.70000000              High     82
## 85   7/8/21 8/12/21      0     0     66 2.20000000              High     66
## 86   7/8/21 8/12/21      0     0     51 1.70000000              High     51
## 87   7/8/21 8/12/21      0     1     36 1.20000000              High     37
## 88   7/8/21 8/12/21      1     0     61 2.03333333              High     62
## 89   7/8/21 8/12/21      1     0     17 0.56666667              High     18
## 90   7/8/21 8/12/21      2     1     66 2.20000000              High     69
## 91   7/8/21 8/12/21      0     1    101 3.36666667              High    102
## 92   7/8/21 8/12/21      3     2     97 3.23333333              High    102
## 93   7/8/21 8/12/21      0     8     87 2.90000000              High     95
## 94   7/8/21 8/12/21      0     0    113 3.76666667              High    113
## 95   7/8/21 8/12/21      0     4     46 1.53333333              High     50
## 96   7/8/21 8/12/21      0     0     82 2.73333333              High     82
## 97   7/8/21 8/12/21      0     0     23 0.76666667               Low     23
## 98   7/8/21 8/12/21      1     0     44 1.46666667               Low     45
## 99   7/8/21 8/12/21      0     0     42 1.40000000               Low     42
## 100  7/8/21 8/12/21      0     0     49 1.63333333               Low     49
## 101  7/8/21 8/12/21      1     1     57 1.90000000               Low     59
## 102  7/8/21 8/12/21      0     2     56 1.86666667               Low     58
## 103  7/8/21 8/12/21      0     0     17 0.56666667               Low     17
## 104  7/8/21 8/12/21      1     0     57 1.90000000               Low     58
## 105  7/8/21 8/12/21      2     1    117 3.90000000               Low    120
## 106 7/15/21 8/19/21      6     1     77 2.56666667               Low     84
## 107 7/15/21 8/19/21      0     1     96 3.20000000               Low     97
## 108 7/15/21 8/19/21      2     2     91 3.03333333               Low     95
## 109 7/15/21 8/19/21      2     1    102 3.40000000               Low    105
## 110 7/15/21 8/19/21      1     3     72 2.40000000               Low     76
## 111 7/15/21 8/19/21      3     1     51 1.70000000               Low     55
## 112 7/15/21 8/19/21      0     0      0 0.00000000               Low      0
## 113 7/15/21 8/19/21      3     0     45 1.50000000               Low     48
## 114 7/15/21 8/19/21      0     5     46 1.53333333               Low     51
## 115 7/15/21 8/19/21      0     0     13 0.43333333               Low     13
## 116 7/15/21 8/19/21      0     1     60 2.00000000               Low     61
## 117 7/15/21 8/19/21      0     1     32 1.06666667               Low     33
## 118 7/15/21 8/19/21      0     2     71 2.36666667               Low     73
## 119 7/15/21 8/19/21      0     1     68 2.26666667               Low     69
## 120  7/8/21 8/12/21      0     1     51 1.70000000               Low     52
## 121  7/8/21 8/12/21      1     0    102 3.40000000               Low    103
## 122  7/8/21 8/12/21      0     0     16 0.53333333               Low     16
## 123  7/8/21 8/12/21      0     0     42 1.40000000               Low     42
## 124 7/15/21 8/19/21      0     0     59 1.96666667               Low     59
## 125 7/15/21 8/19/21      0     0     18 0.60000000               Low     18
## 126 7/15/21 8/19/21      2     0     98 3.26666667               Low    100
## 127 7/15/21 8/19/21      0     0     54 1.80000000               Low     54
## 128 7/15/21 8/19/21      5     1    102 3.40000000               Low    108
## 129 7/15/21 8/19/21      0     2     62 2.06666667               Low     64
## 130 7/22/21 8/26/21      1     3     92 3.06666667               Low     96
## 131 7/22/21 8/26/21      0     0      5 0.16666667               Low      5
## 132  7/8/21 8/12/21      0     2     72 2.40000000               Low     74
## 133  7/8/21 8/12/21      1     1     76 2.53333333               Low     78
## 134  7/8/21 8/12/21      1     0     40 1.33333333               Low     41
## 135 7/22/21 8/26/21      0     0      0 0.00000000               Low      0
## 136 7/22/21 8/26/21      2     0     39 1.30000000               Low     41
## 137 7/22/21 8/26/21      0     1     42 1.40000000               Low     43
## 138 7/22/21 8/26/21      0     0     38 1.26666667               Low     38
## 139 7/22/21 8/26/21      2     0     34 1.13333333               Low     36
## 140 7/22/21 8/26/21      4     1     75 2.50000000               Low     80
## 141  7/8/21 8/12/21      1     0     48 1.60000000               Low     49
## 142  7/8/21 8/12/21      2     2     40 1.33333333               Low     44
## 143  7/8/21 8/12/21      3     4     55 1.83333333               Low     62
## 144  7/8/21 8/12/21      1     5     36 1.20000000               Low     42
## 145  7/8/21 8/12/21      1     0     14 0.46666667               Low     15
## 146  7/8/21 8/12/21      0     0     55 1.83333333               Low     55
## 147  7/8/21 8/12/21      0     0     28 0.93333333               Low     28
## 148  7/8/21 8/12/21      1     2     55 1.83333333               Low     58
## 149  7/8/21 8/12/21      0     0     43 1.43333333               Low     43
## 150  7/8/21 8/12/21      2     1     45 1.50000000               Low     48
## 151  7/8/21 8/12/21     39     6     32 1.06666667               Low     77
## 152  7/8/21 8/12/21      1     2     65 2.16666667               Low     68
## 153  7/8/21 8/12/21      1     0     50 1.66666667               Low     51
## 154  7/8/21 8/12/21      1     1     28 0.93333333               Low     30
## 155  7/8/21 8/12/21      1     0     51 1.70000000               Low     52
## 156  7/8/21 8/12/21      1     0     55 1.83333333               Low     56
## 157  7/8/21 8/12/21      6     1     32 1.06666667            Medium     39
## 158 7/15/21 8/19/21      0     2     43 1.43333333            Medium     45
## 159 7/15/21 8/19/21      0     0     63 2.10000000            Medium     63
## 160 7/15/21 8/19/21      3     0     78 2.60000000            Medium     81
## 161 7/15/21 8/19/21      0     0     60 2.00000000            Medium     60
## 162 7/15/21 8/19/21      0     0     61 2.03333333            Medium     61
## 163 7/15/21 8/19/21      0     0     35 1.16666667            Medium     35
## 164 7/15/21 8/19/21      0     1    106 3.53333333            Medium    107
## 165 7/15/21 8/19/21      3     1     83 2.76666667            Medium     87
## 166 7/15/21 8/19/21      1     1     92 3.06666667            Medium     94
## 167 7/15/21 8/19/21      1     1     26 0.86666667            Medium     28
## 168 7/15/21 8/19/21      1     1     54 1.80000000            Medium     56
## 169 7/15/21 8/19/21      0     0     20 0.66666667            Medium     20
## 170 7/15/21 8/19/21      0     0     41 1.36666667            Medium     41
## 171 7/15/21 8/19/21      1     0     16 0.53333333            Medium     17
## 172 7/22/21 8/26/21      0     0      9 0.30000000            Medium      9
## 173 7/22/21 8/26/21      0     0      0 0.00000000            Medium      0
## 174 7/22/21 8/26/21      0     0      9 0.30000000            Medium      9
## 175 7/22/21 8/26/21      0     0      0 0.00000000            Medium      0
## 176 7/22/21 8/26/21      0     1      5 0.16666667            Medium      6
## 177  7/8/21 8/12/21      0     3     77 2.56666667            Medium     80
## 178  7/8/21 8/12/21      0     0     76 2.53333333            Medium     76
## 179  7/8/21 8/12/21      0     0     30 1.00000000            Medium     30
## 180 7/22/21 8/26/21     18     4     66 2.20000000            Medium     88
## 181 7/22/21 8/26/21      2     0     58 1.93333333            Medium     60
## 182 7/22/21 8/26/21      2     0     98 3.26666667            Medium    100
## 183 7/22/21 8/26/21      0     0     37 1.23333333            Medium     37
## 184  7/8/21 8/12/21      6     0     80 2.66666667            Medium     86
## 185  7/8/21 8/12/21      0     0     42 1.40000000            Medium     42
## 186  7/8/21 8/12/21      0     5     31 1.03333333            Medium     36
## 187  7/8/21 8/12/21      0     2     36 1.20000000            Medium     38
## 188  7/8/21 8/12/21      0     0      1 0.03333333            Medium      1
## 189  7/8/21 8/12/21      0     0     23 0.76666667            Medium     23
## 190  7/8/21 8/12/21      0     2     64 2.13333333            Medium     66
## 191  7/8/21 8/12/21      1     2     32 1.06666667            Medium     35
## 192  7/8/21 8/12/21      1     2     47 1.56666667            Medium     50
## 193  7/8/21 8/12/21      1     2    121 4.03333333            Medium    124
## 194  7/8/21 8/12/21      0     3     79 2.63333333            Medium     82
## 195  7/8/21 8/12/21      1     0     42 1.40000000            Medium     43
## 196  7/8/21 8/12/21      0     2     68 2.26666667            Medium     70
## 197  7/8/21 8/12/21      2     1      8 0.26666667            Medium     11
## 198  7/8/21 8/12/21      1     0     44 1.46666667            Medium     45
## 199 7/15/21 8/19/21      0     2     55 1.83333333            Medium     57
## 200 7/15/21 8/19/21      0     2     63 2.10000000            Medium     65
## 201 7/15/21 8/19/21      0     1     45 1.50000000            Medium     46
## 202 7/15/21 8/19/21      0     2     48 1.60000000            Medium     50
## 203    <NA>    <NA>     NA    NA     NA         NA              <NA>     NA
## 204    <NA>    <NA>     NA    NA     NA         NA              <NA>     NA
## 205    <NA>    <NA>     NA    NA     NA         NA              <NA>     NA
## 206    <NA>    <NA>     NA    NA     NA         NA              <NA>     NA
## 207    <NA>    <NA>     NA    NA     NA         NA              <NA>     NA
## 208    <NA>    <NA>     NA    NA     NA         NA              <NA>     NA
## 209    <NA>    <NA>     NA    NA     NA         NA              <NA>     NA
## 210    <NA>    <NA>     NA    NA     NA         NA              <NA>     NA
## 211    <NA>    <NA>     NA    NA     NA         NA              <NA>     NA
## 212    <NA>    <NA>     NA    NA     NA         NA              <NA>     NA
## 213    <NA>    <NA>     NA    NA     NA         NA              <NA>     NA
##      p_adults     p_pupae    p_larvae He_remain He_remain_exp He_remain_end
## 1   0.8899083 0.082568807 0.027522936 0.9986008     0.9719424     0.9705825
## 2   0.8000000 0.100000000 0.100000000 0.9986008     0.9719424     0.9705825
## 3   0.9719626 0.018691589 0.009345794 0.9986008     0.9825008     0.9811261
## 4   0.9655172 0.011494253 0.022988506 0.9986008     0.9825008     0.9811261
## 5   0.9767442 0.000000000 0.023255814 0.9986008     0.9825008     0.9811261
## 6   0.9880952 0.000000000 0.011904762 0.9986008     0.9806907     0.9793186
## 7   0.9600000 0.000000000 0.040000000 0.9986008     0.9806907     0.9793186
## 8   1.0000000 0.000000000 0.000000000 0.9986008     0.9806907     0.9793186
## 9   0.9772727 0.000000000 0.022727273 0.9986008     0.9806907     0.9793186
## 10  0.9545455 0.000000000 0.045454545 0.9986008     0.9386647     0.9373514
## 11  0.8769231 0.000000000 0.123076923 0.9986008     0.9386647     0.9373514
## 12  0.9821429 0.008928571 0.008928571 0.9986008     0.9797579     0.9783870
## 13  0.9782609 0.000000000 0.021739130 0.9986008     0.9797579     0.9783870
## 14  1.0000000 0.000000000 0.000000000 0.9986008     0.9797579     0.9783870
## 15  0.9473684 0.000000000 0.052631579 0.9986008     0.9797579     0.9783870
## 16  0.9491525 0.033898305 0.016949153 0.9986008     0.9797579     0.9783870
## 17  0.9821429 0.017857143 0.000000000 0.9986008     0.9513318     0.9500007
## 18  0.9090909 0.036363636 0.054545455 0.9986008     0.9803784     0.9790066
## 19  1.0000000 0.000000000 0.000000000 0.9986008     0.9803784     0.9790066
## 20  0.9375000 0.062500000 0.000000000 0.9986008     0.9803784     0.9790066
## 21  0.8800000 0.080000000 0.040000000 0.9986008     0.9836255     0.9822492
## 22  0.9893617 0.010638298 0.000000000 0.9986008     0.9836255     0.9822492
## 23  0.9824561 0.008771930 0.008771930 0.9986008     0.9836255     0.9822492
## 24  0.9763314 0.000000000 0.023668639 0.9986008     0.9802809     0.9789093
## 25  0.9864865 0.013513514 0.000000000 0.9986008     0.9802809     0.9789093
## 26  0.9390244 0.000000000 0.060975610 0.9986008     0.9827964     0.9814213
## 27  0.9767442 0.011627907 0.011627907 0.9986008     0.9827964     0.9814213
## 28  1.0000000 0.000000000 0.000000000 0.9986008     0.9827964     0.9814213
## 29  0.9411765 0.019607843 0.039215686 0.9986008     0.9827964     0.9814213
## 30  1.0000000 0.000000000 0.000000000 0.9986008     0.9827964     0.9814213
## 31  1.0000000 0.000000000 0.000000000 0.9986008     0.9747167     0.9733529
## 32  0.9450549 0.043956044 0.010989011 0.9986008     0.9747167     0.9733529
## 33  0.9680000 0.016000000 0.016000000 0.9986008     0.9747167     0.9733529
## 34  0.9574468 0.028368794 0.014184397 0.9986008     0.9739310     0.9725683
## 35  1.0000000 0.000000000 0.000000000 0.9986008     0.9739310     0.9725683
## 36  0.9838710 0.000000000 0.016129032 0.9986008     0.9852490     0.9838705
## 37  1.0000000 0.000000000 0.000000000 0.9986008     0.9852490     0.9838705
## 38  1.0000000 0.000000000 0.000000000 0.9986008     0.9852490     0.9838705
## 39  1.0000000 0.000000000 0.000000000 0.9986008     0.9852490     0.9838705
## 40  0.9263158 0.052631579 0.021052632 0.9986008     0.9822146     0.9808403
## 41  0.8750000 0.025000000 0.100000000 0.9986008     0.9822146     0.9808403
## 42  0.9702970 0.009900990 0.019801980 0.9986008     0.9855751     0.9841961
## 43  0.9220779 0.064935065 0.012987013 0.9986008     0.9855751     0.9841961
## 44  0.9516129 0.016129032 0.032258065 0.9986008     0.9855751     0.9841961
## 45  0.9347826 0.021739130 0.043478261 0.9986008     0.9855751     0.9841961
## 46  0.9891304 0.000000000 0.010869565 0.9986008     0.9855751     0.9841961
## 47  1.0000000 0.000000000 0.000000000 0.9986008     0.9287849     0.9274854
## 48  0.9102564 0.038461538 0.051282051 0.9986008     0.9751338     0.9737695
## 49  0.6403509 0.035087719 0.324561404 0.9986008     0.9751338     0.9737695
## 50  0.9605263 0.026315789 0.013157895 0.9986008     0.9751338     0.9737695
## 51  0.9021739 0.076086957 0.021739130 0.9986008     0.9751338     0.9737695
## 52  1.0000000 0.000000000 0.000000000 0.9986008     0.9829430     0.9815677
## 53  0.9809524 0.009523810 0.009523810 0.9986008     0.9829430     0.9815677
## 54  0.9420290 0.043478261 0.014492754 0.9986008     0.9829430     0.9815677
## 55  0.9868421 0.000000000 0.013157895 0.9986008     0.9829430     0.9815677
## 56  0.9923077 0.000000000 0.007692308 0.9986008     0.9829430     0.9815677
## 57  1.0000000 0.000000000 0.000000000 0.9986008     0.9829430     0.9815677
## 58  0.9920000 0.000000000 0.008000000 0.9986008     0.9794514     0.9780809
## 59  0.8181818 0.000000000 0.181818182 0.9986008     0.9794514     0.9780809
## 60  1.0000000 0.000000000 0.000000000 0.9986008     0.9794514     0.9780809
## 61  0.9645390 0.014184397 0.021276596 0.9986008     0.9158403     0.9145589
## 62  0.9862069 0.000000000 0.013793103 0.9986008     0.9158403     0.9145589
## 63  0.9531250 0.031250000 0.015625000 0.9986008     0.9428291     0.9415099
## 64  0.9473684 0.000000000 0.052631579 0.9986008     0.9428291     0.9415099
## 65  1.0000000 0.000000000 0.000000000 0.9986008     0.9718326     0.9704729
## 66  1.0000000 0.000000000 0.000000000 0.9986008     0.9718326     0.9704729
## 67  0.9859155 0.014084507 0.000000000 0.9986008     0.9718326     0.9704729
## 68  0.9805825 0.009708738 0.009708738 0.9986008     0.9718326     0.9704729
## 69  1.0000000 0.000000000 0.000000000 0.9986008     0.9718326     0.9704729
## 70  0.9901961 0.000000000 0.009803922 0.9986008     0.9809878     0.9796153
## 71  1.0000000 0.000000000 0.000000000 0.9986008     0.9809878     0.9796153
## 72  1.0000000 0.000000000 0.000000000 0.9986008     0.9809878     0.9796153
## 73  0.9868421 0.000000000 0.013157895 0.9986008     0.9809878     0.9796153
## 74  1.0000000 0.000000000 0.000000000 0.9986008     0.9809878     0.9796153
## 75  1.0000000 0.000000000 0.000000000 0.9986008     0.9809878     0.9796153
## 76  1.0000000 0.000000000 0.000000000 0.9986008     0.9835482     0.9821721
## 77  1.0000000 0.000000000 0.000000000 0.9986008     0.9835482     0.9821721
## 78  1.0000000 0.000000000 0.000000000 0.9986008     0.9835482     0.9821721
## 79  1.0000000 0.000000000 0.000000000 0.9986008     0.9835482     0.9821721
## 80  1.0000000 0.000000000 0.000000000 0.9986008     0.9835482     0.9821721
## 81  0.9777778 0.022222222 0.000000000 0.9986008     0.9855894     0.9842104
## 82  1.0000000 0.000000000 0.000000000 0.9986008     0.9855894     0.9842104
## 83  1.0000000 0.000000000 0.000000000 0.9986008     0.9855894     0.9842104
## 84  0.9878049 0.000000000 0.012195122 0.9986008     0.9855894     0.9842104
## 85  1.0000000 0.000000000 0.000000000 0.9986008     0.9855894     0.9842104
## 86  1.0000000 0.000000000 0.000000000 0.9986008     0.9855894     0.9842104
## 87  0.9729730 0.027027027 0.000000000 0.9986008     0.9855894     0.9842104
## 88  0.9838710 0.000000000 0.016129032 0.9986008     0.9717830     0.9704233
## 89  0.9444444 0.000000000 0.055555556 0.9986008     0.9717830     0.9704233
## 90  0.9565217 0.014492754 0.028985507 0.9986008     0.9717830     0.9704233
## 91  0.9901961 0.009803922 0.000000000 0.9986008     0.9682092     0.9668545
## 92  0.9509804 0.019607843 0.029411765 0.9986008     0.9682092     0.9668545
## 93  0.9157895 0.084210526 0.000000000 0.9986008     0.9682092     0.9668545
## 94  1.0000000 0.000000000 0.000000000 0.9986008     0.9682092     0.9668545
## 95  0.9200000 0.080000000 0.000000000 0.9986008     0.9682092     0.9668545
## 96  1.0000000 0.000000000 0.000000000 0.9986008     0.9682092     0.9668545
## 97  1.0000000 0.000000000 0.000000000 0.5456041     0.9801761     0.5347881
## 98  0.9777778 0.000000000 0.022222222 0.5456041     0.9801761     0.5347881
## 99  1.0000000 0.000000000 0.000000000 0.5456041     0.9801761     0.5347881
## 100 1.0000000 0.000000000 0.000000000 0.5456041     0.9801761     0.5347881
## 101 0.9661017 0.016949153 0.016949153 0.5456041     0.9801761     0.5347881
## 102 0.9655172 0.034482759 0.000000000 0.5411689     0.9684567     0.5240986
## 103 1.0000000 0.000000000 0.000000000 0.5375090     0.9279284     0.4987699
## 104 0.9827586 0.000000000 0.017241379 0.5487325     0.9725983     0.5336963
## 105 0.9750000 0.008333333 0.016666667 0.5487325     0.9725983     0.5336963
## 106 0.9166667 0.011904762 0.071428571 0.5546865     0.9630618     0.5341974
## 107 0.9896907 0.010309278 0.000000000 0.5507359     0.9403195     0.5178677
## 108 0.9578947 0.021052632 0.021052632 0.5507359     0.9403195     0.5178677
## 109 0.9714286 0.009523810 0.019047619 0.5507359     0.9403195     0.5178677
## 110 0.9473684 0.039473684 0.013157895 0.5532470     0.9656184     0.5342255
## 111 0.9272727 0.018181818 0.054545455 0.5532470     0.9656184     0.5342255
## 112       NaN         NaN         NaN 0.5544299     0.6449276     0.3575672
## 113 0.9375000 0.000000000 0.062500000 0.5071881     0.9790690     0.4965722
## 114 0.9019608 0.098039216 0.000000000 0.5071881     0.9790690     0.4965722
## 115 1.0000000 0.000000000 0.000000000 0.4987816     0.9782357     0.4879260
## 116 0.9836066 0.016393443 0.000000000 0.4987816     0.9782357     0.4879260
## 117 0.9696970 0.030303030 0.000000000 0.4987816     0.9782357     0.4879260
## 118 0.9726027 0.027397260 0.000000000 0.4987816     0.9782357     0.4879260
## 119 0.9855072 0.014492754 0.000000000 0.5541380     0.9764292     0.5410765
## 120 0.9807692 0.019230769 0.000000000 0.5448947     0.9741394     0.5308034
## 121 0.9902913 0.000000000 0.009708738 0.5448947     0.9741394     0.5308034
## 122 1.0000000 0.000000000 0.000000000 0.5448947     0.9741394     0.5308034
## 123 1.0000000 0.000000000 0.000000000 0.5448947     0.9741394     0.5308034
## 124 1.0000000 0.000000000 0.000000000 0.5284581     0.9093326     0.4805442
## 125 1.0000000 0.000000000 0.000000000 0.5529282     0.8511338     0.4706159
## 126 0.9800000 0.000000000 0.020000000 0.5435302     0.9711223     0.5278343
## 127 1.0000000 0.000000000 0.000000000 0.5476417     0.9565572     0.5238506
## 128 0.9444444 0.009259259 0.046296296 0.5430646     0.9758713     0.5299612
## 129 0.9687500 0.031250000 0.000000000 0.5430646     0.9758713     0.5299612
## 130 0.9583333 0.031250000 0.010416667 0.5509294     0.9167235     0.5050499
## 131 1.0000000 0.000000000 0.000000000 0.5431402     0.7694643     0.4179270
## 132 0.9729730 0.027027027 0.000000000 0.5460067     0.9556060     0.5217673
## 133 0.9743590 0.012820513 0.012820513 0.5460067     0.9556060     0.5217673
## 134 0.9756098 0.000000000 0.024390244 0.5460067     0.9556060     0.5217673
## 135       NaN         NaN         NaN 0.5523333     0.6083089     0.3359893
## 136 0.9512195 0.000000000 0.048780488 0.5463114     0.9631949     0.5262044
## 137 0.9767442 0.023255814 0.000000000 0.5463114     0.9631949     0.5262044
## 138 1.0000000 0.000000000 0.000000000 0.5463114     0.9631949     0.5262044
## 139 0.9444444 0.000000000 0.055555556 0.5463114     0.9631949     0.5262044
## 140 0.9375000 0.012500000 0.050000000 0.5542119     0.8901046     0.4933065
## 141 0.9795918 0.000000000 0.020408163 0.5489205     0.9707940     0.5328887
## 142 0.9090909 0.045454545 0.045454545 0.5489205     0.9707940     0.5328887
## 143 0.8870968 0.064516129 0.048387097 0.5489205     0.9707940     0.5328887
## 144 0.8571429 0.119047619 0.023809524 0.5489205     0.9707940     0.5328887
## 145 0.9333333 0.000000000 0.066666667 0.5489205     0.9707940     0.5328887
## 146 1.0000000 0.000000000 0.000000000 0.5419774     0.9678511     0.5245534
## 147 1.0000000 0.000000000 0.000000000 0.5419774     0.9678511     0.5245534
## 148 0.9482759 0.034482759 0.017241379 0.5419774     0.9678511     0.5245534
## 149 1.0000000 0.000000000 0.000000000 0.5541201     0.9734332     0.5393989
## 150 0.9375000 0.020833333 0.041666667 0.5541201     0.9734332     0.5393989
## 151 0.4155844 0.077922078 0.506493506 0.5541201     0.9734332     0.5393989
## 152 0.9558824 0.029411765 0.014705882 0.5541201     0.9734332     0.5393989
## 153 0.9803922 0.000000000 0.019607843 0.5274650     0.9808799     0.5173798
## 154 0.9333333 0.033333333 0.033333333 0.5274650     0.9808799     0.5173798
## 155 0.9807692 0.000000000 0.019230769 0.5274650     0.9808799     0.5173798
## 156 0.9821429 0.000000000 0.017857143 0.5274650     0.9808799     0.5173798
## 157 0.8205128 0.025641026 0.153846154 0.6596334     0.9422801     0.6215595
## 158 0.9555556 0.044444444 0.000000000 0.6754272     0.9327651     0.6300149
## 159 1.0000000 0.000000000 0.000000000 0.6803810     0.9756391     0.6638064
## 160 0.9629630 0.000000000 0.037037037 0.6803810     0.9756391     0.6638064
## 161 1.0000000 0.000000000 0.000000000 0.6803810     0.9756391     0.6638064
## 162 1.0000000 0.000000000 0.000000000 0.6803810     0.9756391     0.6638064
## 163 1.0000000 0.000000000 0.000000000 0.6831935     0.8504137     0.5809971
## 164 0.9906542 0.009345794 0.000000000 0.6807728     0.9721768     0.6618315
## 165 0.9540230 0.011494253 0.034482759 0.6807728     0.9721768     0.6618315
## 166 0.9787234 0.010638298 0.010638298 0.6807728     0.9721768     0.6618315
## 167 0.9285714 0.035714286 0.035714286 0.6620374     0.9688824     0.6414364
## 168 0.9642857 0.017857143 0.017857143 0.6620374     0.9688824     0.6414364
## 169 1.0000000 0.000000000 0.000000000 0.6620374     0.9688824     0.6414364
## 170 1.0000000 0.000000000 0.000000000 0.6620374     0.9688824     0.6414364
## 171 0.9411765 0.000000000 0.058823529 0.6620374     0.9688824     0.6414364
## 172 1.0000000 0.000000000 0.000000000 0.6827942     0.9861527     0.6733394
## 173       NaN         NaN         NaN 0.6827942     0.9861527     0.6733394
## 174 1.0000000 0.000000000 0.000000000 0.6827942     0.9861527     0.6733394
## 175       NaN         NaN         NaN 0.6827942     0.9861527     0.6733394
## 176 0.8333333 0.166666667 0.000000000 0.6827298     0.8930853     0.6097360
## 177 0.9625000 0.037500000 0.000000000 0.6799082     0.9736031     0.6619607
## 178 1.0000000 0.000000000 0.000000000 0.6799082     0.9736031     0.6619607
## 179 1.0000000 0.000000000 0.000000000 0.6799082     0.9736031     0.6619607
## 180 0.7500000 0.045454545 0.204545455 0.6804075     0.9418041     0.6408105
## 181 0.9666667 0.000000000 0.033333333 0.6768800     0.9677851     0.6550744
## 182 0.9800000 0.000000000 0.020000000 0.6802035     0.9422407     0.6409154
## 183 1.0000000 0.000000000 0.000000000 0.6802035     0.9422407     0.6409154
## 184 0.9302326 0.000000000 0.069767442 0.6819561     0.9777829     0.6668050
## 185 1.0000000 0.000000000 0.000000000 0.6819561     0.9777829     0.6668050
## 186 0.8611111 0.138888889 0.000000000 0.6819561     0.9777829     0.6668050
## 187 0.9473684 0.052631579 0.000000000 0.6826989     0.9622026     0.6568947
## 188 1.0000000 0.000000000 0.000000000 0.6807065     0.4714889     0.3209456
## 189 1.0000000 0.000000000 0.000000000 0.6819044     0.9771934     0.6663525
## 190 0.9696970 0.030303030 0.000000000 0.6819044     0.9771934     0.6663525
## 191 0.9142857 0.057142857 0.028571429 0.6819044     0.9771934     0.6663525
## 192 0.9400000 0.040000000 0.020000000 0.6800863     0.9695908     0.6594054
## 193 0.9758065 0.016129032 0.008064516 0.6800863     0.9695908     0.6594054
## 194 0.9634146 0.036585366 0.000000000 0.6820356     0.9822570     0.6699342
## 195 0.9767442 0.000000000 0.023255814 0.6820356     0.9822570     0.6699342
## 196 0.9714286 0.028571429 0.000000000 0.6820356     0.9822570     0.6699342
## 197 0.7272727 0.090909091 0.181818182 0.6820356     0.9822570     0.6699342
## 198 0.9777778 0.000000000 0.022222222 0.6820356     0.9822570     0.6699342
## 199 0.9649123 0.035087719 0.000000000 0.6804431     0.9819642     0.6681708
## 200 0.9692308 0.030769231 0.000000000 0.6804431     0.9819642     0.6681708
## 201 0.9782609 0.021739130 0.000000000 0.6804431     0.9819642     0.6681708
## 202 0.9600000 0.040000000 0.000000000 0.6804431     0.9819642     0.6681708
## 203        NA          NA          NA 0.5367998     0.8056432     0.4324691
## 204        NA          NA          NA 0.5386585     0.7452523     0.4014365
## 205        NA          NA          NA 0.5540444     0.4950540     0.2742819
## 206        NA          NA          NA 0.5532775     0.7440450     0.4116634
## 207        NA          NA          NA 0.5528880     0.9892872     0.5469651
## 208        NA          NA          NA 0.5427371     0.4855518     0.2635269
## 209        NA          NA          NA 0.5518545     0.9696290     0.5350941
## 210        NA          NA          NA 0.6802331     0.9285014     0.6315974
## 211        NA          NA          NA 0.6825509     0.8944237     0.6104897
## 212        NA          NA          NA 0.6819650     0.4841994     0.3302071
## 213        NA          NA          NA 0.6809168     0.9614965     0.6546991
SUM_POP_He_ALL_extinct<- Rmisc::summarySE(data_phenotyping_withextinct,
                            measurevar = c("Lambda"),
                       groupvars = c("Genetic_Diversity","He_remain", "He_remain_end",
                                     "Population"), 
                       na.rm=TRUE)


#Add extinction variable
SUM_POP_He_ALL_extinct$Extinction_finalstatus <-ifelse(is.na(SUM_POP_He_ALL_extinct$Lambda),
                                        "yes","no") 


#Replace extinct by NA
SUM_POP_He_ALL_extinct$Lambda <- ifelse(is.na(SUM_POP_He_ALL_extinct$Lambda),
                                        0,SUM_POP_He_ALL_extinct$Lambda)

SUM_POP_He_ALL_extinct$Genetic_Diversity <- ifelse(is.na(SUM_POP_He_ALL_extinct$Genetic_Diversity),
                                        SUM_POP_He_ALL_extinct$Population,
                                        SUM_POP_He_ALL_extinct$Genetic_Diversity)

SUM_POP_He_ALL_extinct$Genetic_Diversity[grep("Low",as.character(SUM_POP_He_ALL_extinct$Population))] <- "Low"
SUM_POP_He_ALL_extinct$Genetic_Diversity[grep("Med",as.character(SUM_POP_He_ALL_extinct$Population))] <- "Medium"
SUM_POP_He_ALL_extinct$Genetic_Diversity[grep("High",as.character(SUM_POP_He_ALL_extinct$Population))] <- "High"


#Order levels 
SUM_POP_He_ALL_extinct$Genetic_Diversity <- factor(SUM_POP_He_ALL_extinct$Genetic_Diversity, levels = c("High", "Medium", "Low"))


PLOT_He_withextinct <- ggplot2::ggplot() +
  geom_point(data = SUM_POP_He_ALL_extinct, 
                               aes(x = He_remain_end, y = Lambda, 
                         colour = Genetic_Diversity, 
                         shape = Extinction_finalstatus, 
                         size = N), alpha = 0.7) +
  ylab("Growth rate") +
  xlab("Remaining heterozygosity") +
  scale_shape_manual(values = c(19, 1)) + 
  theme_LO + 
  labs(color="Genetic diversity") +
  labs(shape="Extinct population") +
  labs(size="Replicates") 

PLOT_He_withextinct

# cowplot::save_plot(here::here("figures","Fitness_He_final.pdf"),   PLOT_He_withextinct,
#           base_height = 12/cm(1), base_width = 18/cm(1), dpi = 610)

3 ANALYSIS

3.1 Probability of extinction

############ 
###### Clean dataset
############ 
#Prepare clean dataset
data_proba_extinction <- data_G6[,c("Block","Population","Genetic_Diversity","Extinction_finalstatus")]
data_proba_extinction$Extinction_finalstatus <- as.factor(data_proba_extinction$Extinction_finalstatus)
data_proba_extinction$y <- ifelse(data_proba_extinction$Extinction_finalstatus == "yes", 1, 0)



############ 
###### Analysis
############ 
#Analysis
mod0 <- glm(y ~ Genetic_Diversity + Block,  
      data = data_proba_extinction, family = binomial)


summary(mod0)
## 
## Call:
## glm(formula = y ~ Genetic_Diversity + Block, family = binomial, 
##     data = data_proba_extinction)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.26011  -0.44258  -0.30957  -0.00005   2.47474  
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)   
## (Intercept)              -21.3144  1855.7435  -0.011  0.99084   
## Genetic_DiversityMedium   18.3002  1855.7433   0.010  0.99213   
## Genetic_DiversityLow      18.5062  1855.7432   0.010  0.99204   
## BlockBlock4                0.7402     1.2714   0.582  0.56046   
## BlockBlock5                3.0006     1.1265   2.664  0.00773 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 72.388  on 83  degrees of freedom
## Residual deviance: 46.833  on 79  degrees of freedom
## AIC: 56.833
## 
## Number of Fisher Scoring iterations: 18
#Test genetic diversity effect
mod1 <- glm(y ~ Block ,  
      data = data_proba_extinction, family = binomial)
anova(mod1, mod0, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: y ~ Block
## Model 2: y ~ Genetic_Diversity + Block
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
## 1        81     58.903                        
## 2        79     46.833  2   12.069 0.002394 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#We keep genetic diversity 

# Perform the lrtest 
(A <- logLik(mod1))
## 'log Lik.' -29.45146 (df=3)
(B <- logLik(mod0))
## 'log Lik.' -23.41669 (df=5)
#Since the logLik() function gives more information than the numeric value, use the as.numeric() function to isolate the numeric value
(teststat <- -2 * (as.numeric(A)-as.numeric(B)))
## [1] 12.06954
#df = 5 - 3 = 2
(p.val <- pchisq(teststat, df = 2, lower.tail = FALSE))
## [1] 0.002394043
lmtest::lrtest(mod1, mod0)
## Likelihood ratio test
## 
## Model 1: y ~ Block
## Model 2: y ~ Genetic_Diversity + Block
##   #Df  LogLik Df  Chisq Pr(>Chisq)   
## 1   3 -29.451                        
## 2   5 -23.417  2 12.069   0.002394 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod0 <- glm(y ~ Genetic_Diversity-1 + Block ,  
      data = data_proba_extinction, family = binomial)
summary(mod0)
## 
## Call:
## glm(formula = y ~ Genetic_Diversity - 1 + Block, family = binomial, 
##     data = data_proba_extinction)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.26011  -0.44258  -0.30957  -0.00005   2.47474  
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)   
## Genetic_DiversityHigh    -21.3144  1855.7435  -0.011  0.99084   
## Genetic_DiversityMedium   -3.0142     1.1293  -2.669  0.00760 **
## Genetic_DiversityLow      -2.8082     1.0659  -2.635  0.00843 **
## BlockBlock4                0.7402     1.2714   0.582  0.56046   
## BlockBlock5                3.0006     1.1265   2.664  0.00773 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 116.449  on 84  degrees of freedom
## Residual deviance:  46.833  on 79  degrees of freedom
## AIC: 56.833
## 
## Number of Fisher Scoring iterations: 18
############ 
###### Predicted value
############ 
#Extract probability of extinction 
(p_high <- plogis(-21.3144))
## [1] 5.536989e-10
(p_med <- plogis(-3.0142))
## [1] 0.04678847
(p_low <- plogis(-2.8082))
## [1] 0.05688267
data_predict_extinct <- data.frame(Genetic_Diversity = levels(data_proba_extinction$Genetic_Diversity),
                                   Block = "Block4")

data_predict_extinct$predict <- predict(mod0, type = "response",
        re.form = NA,
        newdata = data_predict_extinct)
        
data_predict_extinct
##   Genetic_Diversity  Block      predict
## 1              High Block4 1.160723e-09
## 2            Medium Block4 9.329490e-02
## 3               Low Block4 1.122446e-01
p_high
## [1] 5.536989e-10
p_med
## [1] 0.04678847
p_low
## [1] 0.05688267
############ 
###### Posthoc
############ 
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean       SE  df asymp.LCL asymp.UCL
##  High              -20.07 1855.743 Inf  -4451.10  4410.961
##  Medium             -1.77    0.650 Inf     -3.32    -0.216
##  Low                -1.56    0.532 Inf     -2.83    -0.290
## 
## Results are averaged over the levels of: Block 
## Results are given on the logit (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate       SE  df z.ratio p.value
##  High - Medium  -18.300 1855.743 Inf -0.010  0.9999 
##  High - Low     -18.506 1855.743 Inf -0.010  0.9999 
##  Medium - Low    -0.206    0.751 Inf -0.274  0.9594 
## 
## Results are averaged over the levels of: Block 
## Results are given on the log odds ratio (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates

3.2 Population size G6

############ 
###### Analysis
############ 
#Analysis
# mod0 <- lme4::glmer(Nb_adults ~ Genetic_Diversity + Block + (1|Population),  
#       data = data_G6, family = "poisson")


mod0 <- lm(log(Nb_adults) ~ Genetic_Diversity + Block, 
             data = data_G6[data_G6$Extinction_finalstatus == "no",])


mod1 <- lm(log(Nb_adults) ~ Block , 
             data = data_G6[data_G6$Extinction_finalstatus == "no",])

anova(mod0, mod1) # Effect of genetic diversity 
## Analysis of Variance Table
## 
## Model 1: log(Nb_adults) ~ Genetic_Diversity + Block
## Model 2: log(Nb_adults) ~ Block
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1     66 29.706                                  
## 2     68 42.118 -2   -12.412 13.789 9.914e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod0 <- lm(log(Nb_adults) ~ Genetic_Diversity + Block, 
             data = data_G6[data_G6$Extinction_finalstatus == "no",])


summary(mod0)
## 
## Call:
## lm(formula = log(Nb_adults) ~ Genetic_Diversity + Block, data = data_G6[data_G6$Extinction_finalstatus == 
##     "no", ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2417 -0.2947  0.1088  0.4381  1.2624 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              5.11532    0.17980  28.451  < 2e-16 ***
## Genetic_DiversityMedium -0.94813    0.20540  -4.616 1.86e-05 ***
## Genetic_DiversityLow    -0.80532    0.18719  -4.302 5.72e-05 ***
## BlockBlock4             -0.02077    0.18447  -0.113   0.9107    
## BlockBlock5             -0.53922    0.21414  -2.518   0.0142 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6709 on 66 degrees of freedom
## Multiple R-squared:  0.3362, Adjusted R-squared:  0.2959 
## F-statistic: 8.356 on 4 and 66 DF,  p-value: 1.623e-05
############ 
###### Posthoc
############ 
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean    SE df lower.CL upper.CL
##  High                4.93 0.130 66     4.61     5.25
##  Medium              3.98 0.159 66     3.59     4.37
##  Low                 4.12 0.136 66     3.79     4.46
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE df t.ratio p.value
##  High - Medium    0.948 0.205 66  4.616  0.0001 
##  High - Low       0.805 0.187 66  4.302  0.0002 
##  Medium - Low    -0.143 0.207 66 -0.690  0.7704 
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates

3.3 Population size over time

####################### 
####################### MODEL INCLUDING GENETIC DIVERSITY WITHOUT GEN 1  
####################### 
###################### REMOVE GEN 1 ##########################33

#If we dont consider Gen=1 
PLOT_Pop_size_average

fit <- lm(log(Nb_adults+1) ~ Generation_Eggs, data = data[data$Generation_Eggs>=2,])
#second degree
fit2 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,2,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#third degree
fit3 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,3,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#fourth degree
fit4 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,4,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#fourth degree
fit5 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,5,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#exp
fit6 <- lm(log(Nb_adults+1)~exp(Generation_Eggs), data = data[data$Generation_Eggs>=2,])
#log
fit7 <- lm(log(Nb_adults+1)~log(Generation_Eggs), data = data[data$Generation_Eggs>=2,])

fit <- glm(Nb_adults ~ Generation_Eggs, 
           data = data[data$Generation_Eggs>=2,], family = "poisson")
fit2 <- glm(Nb_adults~stats::poly(Generation_Eggs,2,raw=TRUE), 
            data = data[data$Generation_Eggs>=2,], family = "poisson")
fit3 <- glm(Nb_adults~stats::poly(Generation_Eggs,3,raw=TRUE), 
            data = data[data$Generation_Eggs>=2,], family = "poisson")
fit4 <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE), 
            data = data[data$Generation_Eggs>=2,], family = "poisson")
fit5 <- glm(Nb_adults~stats::poly(Generation_Eggs,5,raw=TRUE), 
            data = data[data$Generation_Eggs>=2,], family = "poisson")
fit6 <- glm(Nb_adults~exp(Generation_Eggs), 
            data = data[data$Generation_Eggs>=2,], family = "poisson")
fit7 <- glm(Nb_adults~log(Generation_Eggs), 
            data = data[data$Generation_Eggs>=2,], family = "poisson")





#generate range of 50 numbers starting from 30 and ending at 160
# xx <- seq(1,6, length=100)
# plot(data[data$Generation_Eggs>=2,]$Generation_Eggs,
#      log(data[data$Generation_Eggs>=2,]$Nb_adults+1),pch=19,ylim=c(0,8))
# lines(xx, predict(fit, data.frame(Generation_Eggs=xx)), col="red")
# lines(xx, predict(fit2, data.frame(Generation_Eggs=xx)), col="green")
# lines(xx, predict(fit3, data.frame(Generation_Eggs=xx)), col="blue")
# lines(xx, predict(fit4, data.frame(Generation_Eggs=xx)), col="purple")
# lines(xx, predict(fit5, data.frame(Generation_Eggs=xx)), col="yellow")
# lines(xx, predict(fit6, data.frame(Generation_Eggs=xx)), col="orange")
# lines(xx, predict(fit7, data.frame(Generation_Eggs=xx)), col="pink")
xx <- seq(1,6, length=100)
plot(data[data$Generation_Eggs>=2,]$Generation_Eggs,
     data[data$Generation_Eggs>=2,]$Nb_adults,pch=19)
lines(xx, predict(fit, data.frame(Generation_Eggs=xx),type = "response"), col="red")
lines(xx, predict(fit2, data.frame(Generation_Eggs=xx),type = "response"), col="green")
lines(xx, predict(fit3, data.frame(Generation_Eggs=xx),type = "response"), col="blue")
lines(xx, predict(fit4, data.frame(Generation_Eggs=xx),type = "response"), col="purple")
lines(xx, predict(fit5, data.frame(Generation_Eggs=xx),type = "response"), col="yellow")
lines(xx, predict(fit6, data.frame(Generation_Eggs=xx),type = "response"), col="orange")
lines(xx, predict(fit7, data.frame(Generation_Eggs=xx),type = "response"), col="pink")

# Check the r squared
rsq::rsq(fit,adj=TRUE)
## [1] 0.4442203
rsq::rsq(fit2,adj=TRUE)
## [1] 0.5938535
rsq::rsq(fit3,adj=TRUE)
## [1] 0.5928435
rsq::rsq(fit4,adj=TRUE)
## [1] 0.592052
rsq::rsq(fit5,adj=TRUE)
## [1] 0.592052
rsq::rsq(fit6,adj=TRUE)
## [1] 0.121301
rsq::rsq(fit7,adj=TRUE)
## [1] 0.5205104
# Best model is fit 2


data$gen_square <- data$Generation_Eggs^2
fit2 <- glm(Nb_adults~stats::poly(Generation_Eggs,2,raw=TRUE), 
            data = data[data$Generation_Eggs>=2,], family = "poisson")
fit2 <- glm(Nb_adults~Generation_Eggs + gen_square, 
            data = data[data$Generation_Eggs>=2,], family = "poisson")

# fit2 <- lm(log(Nb_adults+1)~poly(Generation_Eggs,2,raw=TRUE), data = data[data$Generation_Eggs>=2,])
# fit2 <- lm(log(Nb_adults+1)~ Generation_Eggs + gen_square, data = data[data$Generation_Eggs>=2,])
#Add the interaction with Genetic diversity
fit2 <- glm(Nb_adults~ Generation_Eggs*Genetic_Diversity +
             gen_square*Genetic_Diversity,
           data = data, family = "poisson")


summary(fit2)
## 
## Call:
## glm(formula = Nb_adults ~ Generation_Eggs * Genetic_Diversity + 
##     gen_square * Genetic_Diversity, family = "poisson", data = data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -15.961   -6.480   -2.662    5.254   21.115  
## 
## Coefficients:
##                                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)                              5.070332   0.026463 191.604  < 2e-16
## Generation_Eggs                          0.156559   0.017811   8.790  < 2e-16
## Genetic_DiversityMedium                 -0.164697   0.042046  -3.917 8.96e-05
## Genetic_DiversityLow                     0.024718   0.037632   0.657 0.511290
## gen_square                              -0.036751   0.002566 -14.324  < 2e-16
## Generation_Eggs:Genetic_DiversityMedium  0.080423   0.029501   2.726 0.006409
## Generation_Eggs:Genetic_DiversityLow    -0.096063   0.026425  -3.635 0.000278
## Genetic_DiversityMedium:gen_square      -0.035166   0.004448  -7.905 2.67e-15
## Genetic_DiversityLow:gen_square         -0.009935   0.003956  -2.511 0.012038
##                                            
## (Intercept)                             ***
## Generation_Eggs                         ***
## Genetic_DiversityMedium                 ***
## Genetic_DiversityLow                       
## gen_square                              ***
## Generation_Eggs:Genetic_DiversityMedium ** 
## Generation_Eggs:Genetic_DiversityLow    ***
## Genetic_DiversityMedium:gen_square      ***
## Genetic_DiversityLow:gen_square         *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 39168  on 486  degrees of freedom
## Residual deviance: 30719  on 478  degrees of freedom
##   (17 observations deleted due to missingness)
## AIC: 33747
## 
## Number of Fisher Scoring iterations: 5
####### Add genetic diversity 

#Test oversdispersion
mod0 <- glm(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity + 
              Block , 
            data = data[data$Generation_Eggs>=2,], family = "poisson")
sqrt(deviance(mod0)/(nobs(mod0)-length(coef(mod0))))
## [1] 6.395475
sigma(mod0)
## [1] 6.395475
# ccl: overdispersion


mod1 <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity + 
                Block  + (1|obs), 
            data = data[data$Generation_Eggs>=2,], family = "poisson")
blmeco::dispersion_glmer(mod1)
## [1] 1.076696
sigma(mod1)
## [1] 1
#Convergence issue 
max(abs(with(mod1@optinfo$derivs, solve(Hessian, gradient)))) #should be <0.001
## [1] 0.004551091
# Test effect
mod1 <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity + 
                Block  + (1|obs), 
            data = data[data$Generation_Eggs>=2,], family = "poisson")
mod2 <- lme4::glmer(Nb_adults~Generation_Eggs + gen_square + Genetic_Diversity + 
                Block  + (1|obs), 
            data = data[data$Generation_Eggs>=2,], family = "poisson")
anova(mod1, mod2, test ="Chisq")  
## Data: data[data$Generation_Eggs >= 2, ]
## Models:
## mod2: Nb_adults ~ Generation_Eggs + gen_square + Genetic_Diversity + 
## mod2:     Block + (1 | obs)
## mod1: Nb_adults ~ Generation_Eggs * Genetic_Diversity + gen_square * 
## mod1:     Genetic_Diversity + Block + (1 | obs)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## mod2    8 4669.1 4701.1 -2326.6   4653.1                        
## mod1   12 4660.8 4708.8 -2318.4   4636.8 16.376  4   0.002554 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Predict data
filldata <- data.frame(Genetic_Diversity = rep(levels(data$Genetic_Diversity),100), 
                       Generation_Eggs = rep(seq(2,6, length = 100),each=3),
                       Block = "Block4",
                       gen_square = (rep(seq(2,6, length = 100),each=3))^2,
                       Estimates = NA)

filldata$Estimates <- predict(mod1, newdata=filldata, type = "response", re.form=~0)


#Plot
PLOT_Pop_size <- ggplot2::ggplot(data = data[data$Generation_Eggs>=2,]) + 
  geom_line(data = filldata, aes(x = Generation_Eggs, 
                                y = Estimates, 
                               colour = Genetic_Diversity), size = 1.3) +
  geom_point(data = data, 
                      aes(x = Generation_Eggs, 
                                y = Nb_adults, 
                                colour = Genetic_Diversity),
             position = position_dodge(0.4), size = 1.4, alpha = 0.6) +
  ylab("Population size") +
  labs(color="Genetic diversity") +
  xlab("Generation") +
  theme_LO
PLOT_Pop_size

emmeans::emmeans(mod1, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean     SE  df asymp.LCL asymp.UCL
##  High                4.88 0.0790 Inf      4.69      5.07
##  Medium              4.19 0.0885 Inf      3.98      4.40
##  Low                 4.04 0.0733 Inf      3.87      4.22
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE  df z.ratio p.value
##  High - Medium    0.688 0.119 Inf 5.786   <.0001 
##  High - Low       0.835 0.107 Inf 7.778   <.0001 
##  Medium - Low     0.147 0.115 Inf 1.287   0.4028 
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates
###################################################################
###################################################################
###################################################################
####################  INCLUDING GEN 1 ###########################
PLOT_Pop_size_average

fit <- glm(Nb_adults ~ Generation_Eggs, data = data, family = "poisson")
#second degree
fit2 <- glm(Nb_adults~stats::poly(Generation_Eggs,2,raw=TRUE), data = data, family = "poisson")
#third degree
fit3 <- glm(Nb_adults~stats::poly(Generation_Eggs,3,raw=TRUE), data = data, family = "poisson")
#fourth degree
fit4 <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE), data = data, family = "poisson")
#fourth degree
fit5 <- glm(Nb_adults~stats::poly(Generation_Eggs,5,raw=TRUE), data = data, family = "poisson")
#exp
fit6 <- glm(Nb_adults~exp(Generation_Eggs), data = data, family = "poisson")
#log
fit7 <- glm(Nb_adults~log(Generation_Eggs), data = data, family = "poisson")


#generate range of 50 numbers starting from 30 and ending at 160
xx <- seq(1,6, length=100)
plot(data$Generation_Eggs,
     data$Nb_adults,pch=19)
lines(xx, predict(fit, data.frame(Generation_Eggs=xx),type = "response"), col="red")
lines(xx, predict(fit2, data.frame(Generation_Eggs=xx),type = "response"), col="green")
lines(xx, predict(fit3, data.frame(Generation_Eggs=xx),type = "response"), col="blue")
lines(xx, predict(fit4, data.frame(Generation_Eggs=xx),type = "response"), col="purple")
lines(xx, predict(fit5, data.frame(Generation_Eggs=xx),type = "response"), col="yellow")
lines(xx, predict(fit6, data.frame(Generation_Eggs=xx),type = "response"), col="orange")
lines(xx, predict(fit7, data.frame(Generation_Eggs=xx),type = "response"), col="pink")

# Check the r squared
summary(fit)$adj.r.squared
## NULL
summary(fit2)$adj.r.squared
## NULL
summary(fit3)$adj.r.squared
## NULL
summary(fit4)$adj.r.squared # best r square
## NULL
summary(fit5)$adj.r.squared 
## NULL
summary(fit6)$adj.r.squared 
## NULL
summary(fit7)$adj.r.squared 
## NULL
rsq::rsq(fit,adj=TRUE)
## [1] 0.1067643
rsq::rsq(fit2,adj=TRUE)
## [1] 0.1486733
rsq::rsq(fit3,adj=TRUE)
## [1] 0.5240033
rsq::rsq(fit4,adj=TRUE)
## [1] 0.5980368
rsq::rsq(fit5,adj=TRUE)
## [1] 0.5989426
rsq::rsq(fit6,adj=TRUE)
## [1] 0.07185974
rsq::rsq(fit7,adj=TRUE)
## [1] 0.04937507
#
fit4 <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE), data = data, family = "poisson")



####################### 
####################### MODEL INCLUDING GENETIC DIVERSITY and GEN 1 
####################### 

#Test oversdispersion
fit4 <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity + (1|Block), data = data, family = "poisson")
blmeco::dispersion_glmer(fit4)
## [1] 5.760264
mod <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity + 
               (1|Block) + (1|obs), data = data, family = "poisson")
blmeco::dispersion_glmer(mod)
## [1] 1.083377
#Convergence issue 
max(abs(with(mod@optinfo$derivs, solve(Hessian, gradient)))) #should be <0.001
## [1] 0.06384154
#Predict data
filldata <- data.frame(Genetic_Diversity = rep(levels(data$Genetic_Diversity),100), 
                       Generation_Eggs = rep(seq(1,6, length = 100),each=3),
                       Block = "Block4", 
                       Estimates = NA)
filldata$Estimates <- predict(fit4, newdata=filldata, type = "response")


#Plot
PLOT_Pop_size <- ggplot2::ggplot(data = data) + 
  geom_line(data = filldata, aes(x = Generation_Eggs, 
                                y = Estimates, 
                               colour = Genetic_Diversity), size = 1.3) +
  geom_point(data = data, 
                      aes(x = Generation_Eggs, 
                                y = Nb_adults, 
                                colour = Genetic_Diversity),
             position = position_dodge(0.4), size = 1.4, alpha = 0.6) +
  ylab("Population size") +
  labs(color="Genetic diversity") +
  xlab("Generation") +
  theme_LO
PLOT_Pop_size

# fit2 <- lm(log(Nb_adults+1)~ Generation_Eggs*Genetic_Diversity + 
#              gen_square*Genetic_Diversity, 
#            data = data[data$Generation_Eggs>=2,])
# fit2_test <- lm(log(Nb_adults+1)~ Generation_Eggs+Genetic_Diversity + 
#              gen_square, 
#            data = data[data$Generation_Eggs>=2,])

mod <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity + 
               (1|Block) + (1|obs), data = data, family = "poisson")
mod1 <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)+Genetic_Diversity + 
               (1|Block) + (1|obs), data = data, family = "poisson")
anova(mod,mod1, test = "Chisq")
## Data: data
## Models:
## mod1: Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) + Genetic_Diversity + 
## mod1:     (1 | Block) + (1 | obs)
## mod: Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity + 
## mod:     (1 | Block) + (1 | obs)
##      npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)    
## mod1    9 5611.2 5648.9 -2796.6   5593.2                        
## mod    17 5594.8 5666.0 -2780.4   5560.8 32.47  8  7.672e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Signifcant interaction 
summary(mod)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity +  
##     (1 | Block) + (1 | obs)
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##   5594.8   5666.0  -2780.4   5560.8      470 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.03685 -0.03969  0.00800  0.05695  0.26801 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  obs    (Intercept) 0.64646  0.8040  
##  Block  (Intercept) 0.07945  0.2819  
## Number of obs: 487, groups:  obs, 487; Block, 3
## 
## Fixed effects:
##                                                                       Estimate
## (Intercept)                                                          -1.993477
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         10.858800
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         -5.147603
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                          0.942314
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         -0.058998
## Genetic_DiversityMedium                                              -0.352704
## Genetic_DiversityLow                                                 -0.787854
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium  0.772686
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium -0.503761
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium  0.095237
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium -0.006018
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow     1.489498
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow    -0.833931
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow     0.138807
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow    -0.007028
##                                                                      Std. Error
## (Intercept)                                                            1.238244
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                           1.985855
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                           1.032426
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                           0.213266
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                           0.015171
## Genetic_DiversityMedium                                                1.833248
## Genetic_DiversityLow                                                   1.630486
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium   2.976126
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium   1.553127
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium   0.322073
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium   0.022992
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow      2.637956
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow      1.373155
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow      0.284239
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow      0.020268
##                                                                      z value
## (Intercept)                                                           -1.610
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                           5.468
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                          -4.986
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                           4.418
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                          -3.889
## Genetic_DiversityMedium                                               -0.192
## Genetic_DiversityLow                                                  -0.483
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium   0.260
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium  -0.324
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium   0.296
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium  -0.262
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow      0.565
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow     -0.607
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow      0.488
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow     -0.347
##                                                                      Pr(>|z|)
## (Intercept)                                                          0.107415
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         4.55e-08
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         6.17e-07
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                         9.94e-06
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         0.000101
## Genetic_DiversityMedium                                              0.847434
## Genetic_DiversityLow                                                 0.628952
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 0.795151
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 0.745671
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.767459
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.793538
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow    0.572318
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow    0.543645
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow    0.625306
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow    0.728786
##                                                                         
## (Intercept)                                                             
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         ***
## Genetic_DiversityMedium                                                 
## Genetic_DiversityLow                                                    
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model failed to converge with max|grad| = 0.398806 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
summary(mod1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) + Genetic_Diversity +  
##     (1 | Block) + (1 | obs)
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##   5611.2   5648.9  -2796.6   5593.2      478 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.05397 -0.05185  0.01426  0.07054  0.21832 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  obs    (Intercept) 0.69318  0.8326  
##  Block  (Intercept) 0.07665  0.2769  
## Number of obs: 487, groups:  obs, 487; Block, 3
## 
## Fixed effects:
##                                              Estimate Std. Error z value
## (Intercept)                                  -2.01996    0.81642  -2.474
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 11.74179    1.30992   8.964
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 -5.66290    0.68789  -8.232
## stats::poly(Generation_Eggs, 4, raw = TRUE)3  1.03332    0.14306   7.223
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 -0.06412    0.01022  -6.274
## Genetic_DiversityMedium                      -0.55984    0.10012  -5.592
## Genetic_DiversityLow                         -0.67366    0.09048  -7.445
##                                              Pr(>|z|)    
## (Intercept)                                    0.0134 *  
## stats::poly(Generation_Eggs, 4, raw = TRUE)1  < 2e-16 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)2  < 2e-16 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 5.08e-13 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 3.52e-10 ***
## Genetic_DiversityMedium                      2.25e-08 ***
## Genetic_DiversityLow                         9.69e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                    (Intr) s::(G_E,4,r=TRUE)1 s::(G_E,4,r=TRUE)2
## s::(G_E,4,r=TRUE)1 -0.965                                      
## s::(G_E,4,r=TRUE)2  0.941 -0.993                               
## s::(G_E,4,r=TRUE)3 -0.915  0.976             -0.995            
## s::(G_E,4,r=TRUE)4  0.889 -0.956              0.984            
## Gntc_DvrstM        -0.062  0.007             -0.008            
## Gntc_DvrstL        -0.065  0.004             -0.005            
##                    s::(G_E,4,r=TRUE)3 s::(G_E,4,r=TRUE)4 Gnt_DM
## s::(G_E,4,r=TRUE)1                                             
## s::(G_E,4,r=TRUE)2                                             
## s::(G_E,4,r=TRUE)3                                             
## s::(G_E,4,r=TRUE)4 -0.997                                      
## Gntc_DvrstM         0.008             -0.009                   
## Gntc_DvrstL         0.005             -0.006              0.494
## convergence code: 0
## Model failed to converge with max|grad| = 0.053402 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
emmeans::emmeans(mod, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean    SE  df asymp.LCL asymp.UCL
##  High                4.55 0.206 Inf      4.06      5.04
##  Medium              3.93 0.216 Inf      3.42      4.45
##  Low                 3.69 0.201 Inf      3.21      4.17
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE  df z.ratio p.value
##  High - Medium    0.619 0.187 Inf 3.310   0.0027 
##  High - Low       0.861 0.168 Inf 5.115   <.0001 
##  Medium - Low     0.242 0.184 Inf 1.319   0.3844 
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates
####################### 
####################### MODEL INCLUDING ONLY RESCUED POPULATIONS
####################### 
#Same but only for population not extinct
fit_fin <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity,
            data = data[data$Extinction_finalstatus=="no",], family = "poisson")

filldata$Estimates_Withoutextpop <- predict(fit_fin, newdata=filldata, type = "response")

PLOT_Pop_size <- ggplot2::ggplot(data = data) + 
  geom_line(data = filldata, aes(x = Generation_Eggs, 
                                y = Estimates_Withoutextpop, 
                               colour = Genetic_Diversity), size = 1.5, linetype = "longdash") +
  geom_point(data = data[data$Extinction_finalstatus=="no",], 
                      aes(x = Generation_Eggs, 
                                y = Nb_adults, 
                                colour = Genetic_Diversity),
             position = position_dodge(0.4), size = 1.4, alpha = 0.5) +
  ylab("Population size") +
  labs(color="Genetic diversity") +
  xlab("Generation") +
  theme_LO
PLOT_Pop_size

mod <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity + 
               (1|Block) + (1|obs), 
             data = data[data$Extinction_finalstatus=="no",], family = "poisson")
mod1 <- lme4::glmer(Nb_adults~poly(Generation_Eggs,4,raw=TRUE)+Genetic_Diversity + 
               (1|Block) + (1|obs), 
              data = data[data$Extinction_finalstatus=="no",], family = "poisson")
anova(mod,mod1, test = "Chisq")
## Data: data[data$Extinction_finalstatus == "no", ]
## Models:
## mod1: Nb_adults ~ poly(Generation_Eggs, 4, raw = TRUE) + Genetic_Diversity + 
## mod1:     (1 | Block) + (1 | obs)
## mod: Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity + 
## mod:     (1 | Block) + (1 | obs)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## mod1    9 4759.9 4796.3 -2370.9   4741.9                       
## mod    17 4756.3 4825.1 -2361.1   4722.3 19.599  8    0.01196 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity +  
##     (1 | Block) + (1 | obs)
##    Data: data[data$Extinction_finalstatus == "no", ]
## 
##      AIC      BIC   logLik deviance df.resid 
##   4756.3   4825.1  -2361.1   4722.3      407 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.57063 -0.04879  0.00393  0.06967  0.29650 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  obs    (Intercept) 0.36624  0.6052  
##  Block  (Intercept) 0.02456  0.1567  
## Number of obs: 424, groups:  obs, 424; Block, 3
## 
## Fixed effects:
##                                                                        Estimate
## (Intercept)                                                          -1.9250178
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         10.7535583
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         -5.0922588
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                          0.9316535
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         -0.0583188
## Genetic_DiversityMedium                                               1.1892488
## Genetic_DiversityLow                                                  0.2794991
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium -2.0341925
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium  1.0112012
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium -0.2039827
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium  0.0137179
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow    -0.4034287
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow     0.0777534
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow    -0.0101264
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow     0.0006044
##                                                                      Std. Error
## (Intercept)                                                           0.9705855
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                          1.5722049
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                          0.8204596
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                          0.1698215
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                          0.0120893
## Genetic_DiversityMedium                                               1.5049889
## Genetic_DiversityLow                                                  1.3816159
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium  2.4448787
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium  1.2751857
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium  0.2640143
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium  0.0188078
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow     2.2495506
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow     1.1751948
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow     0.2435505
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow     0.0173582
##                                                                      z value
## (Intercept)                                                           -1.983
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                           6.840
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                          -6.207
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                           5.486
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                          -4.824
## Genetic_DiversityMedium                                                0.790
## Genetic_DiversityLow                                                   0.202
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium  -0.832
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium   0.793
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium  -0.773
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium   0.729
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow     -0.179
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow      0.066
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow     -0.042
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow      0.035
##                                                                      Pr(>|z|)
## (Intercept)                                                            0.0473
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         7.93e-12
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         5.41e-10
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                         4.11e-08
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         1.41e-06
## Genetic_DiversityMedium                                                0.4294
## Genetic_DiversityLow                                                   0.8397
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium   0.4054
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium   0.4278
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium   0.4397
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium   0.4658
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow      0.8577
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow      0.9472
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow      0.9668
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow      0.9722
##                                                                         
## (Intercept)                                                          *  
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         ***
## Genetic_DiversityMedium                                                 
## Genetic_DiversityLow                                                    
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model failed to converge with max|grad| = 3.25364 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
emmeans::emmeans(mod, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean    SE  df asymp.LCL asymp.UCL
##  High                4.53 0.133 Inf      4.21      4.85
##  Medium              4.30 0.151 Inf      3.94      4.66
##  Low                 4.01 0.137 Inf      3.68      4.33
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE  df z.ratio p.value
##  High - Medium    0.230 0.153 Inf 1.507   0.2874 
##  High - Low       0.523 0.140 Inf 3.746   0.0005 
##  Medium - Low     0.293 0.157 Inf 1.861   0.1503 
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates
#Signifcant interaction 

3.4 Growth rate G6

head(data_phenotyping)
##   Population   Week  Block ID_Rep Divsersity Population.1 Box Initial.N   Start
## 1      High1 5-week Block3     52       High            1   1        30  7/8/21
## 2      High1 5-week Block3     53       High            1   2        30  7/8/21
## 3     High13 5-week Block4    130       High           13   2        30 7/15/21
## 4     High13 5-week Block4    131       High           13   3        30 7/15/21
## 5     High13 5-week Block4    129       High           13   1        30 7/15/21
## 6     High14 5-week Block4    132       High           14   1        30 7/15/21
##       End Larvae Pupae Adults    Lambda Genetic_Diversity Nb_ttx  p_adults
## 1 8/12/21      3     9     97 3.2333333              High    109 0.8899083
## 2 8/12/21      1     1      8 0.2666667              High     10 0.8000000
## 3 8/19/21      1     2    104 3.4666667              High    107 0.9719626
## 4 8/19/21      2     1     84 2.8000000              High     87 0.9655172
## 5 8/19/21      2     0     84 2.8000000              High     86 0.9767442
## 6 8/19/21      1     0     83 2.7666667              High     84 0.9880952
##      p_pupae    p_larvae He_remain He_remain_exp He_remain_end
## 1 0.08256881 0.027522936 0.9986008     0.9719424     0.9705825
## 2 0.10000000 0.100000000 0.9986008     0.9719424     0.9705825
## 3 0.01869159 0.009345794 0.9986008     0.9825008     0.9811261
## 4 0.01149425 0.022988506 0.9986008     0.9825008     0.9811261
## 5 0.00000000 0.023255814 0.9986008     0.9825008     0.9811261
## 6 0.00000000 0.011904762 0.9986008     0.9806907     0.9793186
#Without considering extinct populations 

#tapply(data_phenotyping$Lambda,data_phenotyping$Population,length)




#############################################################
################## Determine distribution ###################
#############################################################

## Poisson lognormal model
mlog <- lme4::lmer(log(Lambda+1) ~ Genetic_Diversity + Block + (1|Population), 
                   data=data_phenotyping)
#summary(mlog)

# Square root
msqrt <- lme4::lmer(sqrt(Lambda) ~ Genetic_Diversity + Block + (1|Population), 
                   data=data_phenotyping)
#summary(msqrt)

# Normal
mnorm <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population), 
                   data=data_phenotyping)


## Compare AIC
AIC(mlog, msqrt,mnorm)
##       df      AIC
## mlog   7 169.8536
## msqrt  7 215.0823
## mnorm  7 626.7241
#Square root a better fits to the data 
#But we compare different values: Growth rate and Growth rate+1

# ## Simulate data 
x.teo.log <- unlist(simulate(mlog))
x.teo.sqrt <- unlist(simulate(msqrt))
x.teo.norm <- unlist(simulate(mnorm))

# ## QQplot to compared Negative binomial and Poisson log normal distributions
qqplot(data_phenotyping$Lambda, x.teo.log,
       main="QQ-plot", xlab="Observed data", ylab="Simulated data",
       las=1, xlim = c(0, 10), ylim = c(0, 10), col='blue') ## QQ-plot
points(sort(data_phenotyping$Lambda), sort(x.teo.sqrt), pch=1, col='red')
points(sort(data_phenotyping$Lambda), sort(x.teo.norm), pch=1, col='green')
abline(0,1)
# ## Add legend
legend(0, 7, legend=c("Log", "Sqrt", "Normal"),
       col=c("blue", "red", "green"), 
       pch=1, bty="n")

# ## Normal distribution provides a good fit to the data


# QQ plot
# qqnorm(resid(mlog))
# qqline(resid(mlog))
# qqnorm(resid(msqrt))
# qqline(resid(msqrt))
# qqnorm(resid(mnorm))
# qqline(resid(mnorm)) #best fit

#Distribution Residuals
hist(residuals(mlog),col="yellow",freq=F)

hist(residuals(msqrt),col="yellow",freq=F)  #Seems the better fit

hist(residuals(mnorm),col="yellow",freq=F) #Seems the better fit

#Test normality
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(mnorm), "norm"))
g1$kstest 
##     1-mle-norm 
## "not rejected"
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(msqrt), "norm"))
g1$kstest 
##     1-mle-norm 
## "not rejected"
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(mnorm), "norm"))
g1$kstest
##     1-mle-norm 
## "not rejected"
## Compare AIC
AIC(mlog, msqrt, mnorm)
##       df      AIC
## mlog   7 169.8536
## msqrt  7 215.0823
## mnorm  7 626.7241
#SLog a better fits to the data 





#############################################################
##################        Analysis        ###################
#############################################################

mod0 <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population), 
                   data=data_phenotyping)
summary(mod0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Lambda ~ Genetic_Diversity + Block + (1 | Population)
##    Data: data_phenotyping
## 
## REML criterion at convergence: 612.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4322 -0.6449 -0.0938  0.5243  3.3448 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Population (Intercept) 0.2958   0.5439  
##  Residual               0.7591   0.8713  
## Number of obs: 217, groups:  Population, 77
## 
## Fixed effects:
##                         Estimate Std. Error t value
## (Intercept)               2.5956     0.1904  13.629
## Genetic_DiversityMedium  -1.1229     0.2271  -4.944
## Genetic_DiversityLow     -0.8852     0.2128  -4.159
## BlockBlock4               0.2452     0.2059   1.191
## BlockBlock5              -0.1742     0.2407  -0.724
## 
## Correlation of Fixed Effects:
##             (Intr) Gnt_DM Gnt_DL BlckB4
## Gntc_DvrstM -0.488                     
## Gntc_DvrstL -0.559  0.401              
## BlockBlock4 -0.605  0.044  0.075       
## BlockBlock5 -0.593  0.111  0.185  0.457
mod1 <- lme4::lmer(Lambda ~ Block + (1|Population), 
                   data=data_phenotyping)
anova(mod0,mod1,test="Chisq") # difference
## Data: data_phenotyping
## Models:
## mod1: Lambda ~ Block + (1 | Population)
## mod0: Lambda ~ Genetic_Diversity + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## mod1    5 641.07 657.97 -315.54   631.07                         
## mod0    7 618.35 642.01 -302.18   604.35 26.718  2  1.579e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############ 
###### Posthoc
############ 
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean    SE   df lower.CL upper.CL
##  High                2.62 0.136 52.2     2.28     2.96
##  Medium              1.50 0.183 71.3     1.05     1.94
##  Low                 1.73 0.165 83.1     1.33     2.14
## 
## Results are averaged over the levels of: Block 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE   df t.ratio p.value
##  High - Medium    1.123 0.228 63.5  4.931  <.0001 
##  High - Low       0.885 0.213 69.5  4.147  0.0003 
##  Medium - Low    -0.238 0.242 74.8 -0.983  0.5901 
## 
## Results are averaged over the levels of: Block 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates
############ 
###### Predict
############
data_predict_extinct <- data.frame(Genetic_Diversity = levels(data_phenotyping$Genetic_Diversity),
                                   Block = "Block4")

data_predict_extinct$predict <- predict(mod0, 
                                        type = "response",
                                        re.form = NA,
                                        newdata = data_predict_extinct)

3.5 Remaining He difference

Rmisc::summarySE(data_G2,
          measurevar = c("He_remain"),
          groupvars = c("Genetic_Diversity"),
          na.rm=TRUE)
##   Genetic_Diversity  N He_remain          sd          se          ci
## 1              High 27 0.9986008 0.000000000 0.000000000 0.000000000
## 2            Medium 23 0.6791246 0.006059307 0.001263453 0.002620241
## 3               Low 34 0.5441887 0.012687849 0.002175948 0.004427000
mfull <- lm(He_remain ~ Genetic_Diversity, data=data_G2)
mless <- lm(He_remain ~ 1, data=data_G2)
anova(mfull,mless) #Keep effect
## Analysis of Variance Table
## 
## Model 1: He_remain ~ Genetic_Diversity
## Model 2: He_remain ~ 1
##   Res.Df    RSS Df Sum of Sq     F    Pr(>F)    
## 1     81 0.0061                                 
## 2     83 3.1868 -2   -3.1807 21048 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############ 
###### Posthoc
############ 
emmeans::emmeans(mfull, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean       SE df lower.CL upper.CL
##  High              0.9986 0.001673 81   0.9945   1.0027
##  Medium            0.6791 0.001812 81   0.6747   0.6835
##  Low               0.5442 0.001491 81   0.5406   0.5478
## 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate      SE df t.ratio p.value
##  High - Medium    0.319 0.00247 81 129.527 <.0001 
##  High - Low       0.454 0.00224 81 202.800 <.0001 
##  Medium - Low     0.135 0.00235 81  57.498 <.0001 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates
Rmisc::summarySE(data_G2,
          measurevar = c("He_remain"),
          groupvars = c("Genetic_Diversity"),
          na.rm=TRUE)
##   Genetic_Diversity  N He_remain          sd          se          ci
## 1              High 27 0.9986008 0.000000000 0.000000000 0.000000000
## 2            Medium 23 0.6791246 0.006059307 0.001263453 0.002620241
## 3               Low 34 0.5441887 0.012687849 0.002175948 0.004427000
#Calcul for end:
Rmisc::summarySE(data_G2[data_G2$Extinction_finalstatus!="yes",],
          measurevar = c("He_remain_end"),
          groupvars = c("Genetic_Diversity"),
          na.rm=TRUE)
##   Genetic_Diversity  N He_remain_end         sd          se          ci
## 1              High 27     0.9697741 0.01868272 0.003595491 0.007390637
## 2            Medium 18     0.6482803 0.02453741 0.005783522 0.012202164
## 3               Low 26     0.5146307 0.02759467 0.005411760 0.011145728
mfull <- lm(He_remain_end ~ Genetic_Diversity, data=data_G2[data_G2$Extinction_finalstatus!="yes",])
mless <- lm(He_remain_end ~ 1, data=data_G2[data_G2$Extinction_finalstatus!="yes",])
anova(mfull,mless) #Keep effect
## Analysis of Variance Table
## 
## Model 1: He_remain_end ~ Genetic_Diversity
## Model 2: He_remain_end ~ 1
##   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
## 1     68 0.03835                                  
## 2     70 2.91180 -2   -2.8735 2547.7 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############ 
###### Posthoc
############ 
emmeans::emmeans(mfull, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean      SE df lower.CL upper.CL
##  High               0.970 0.00457 68    0.959    0.981
##  Medium             0.648 0.00560 68    0.635    0.662
##  Low                0.515 0.00466 68    0.503    0.526
## 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate      SE df t.ratio p.value
##  High - Medium    0.321 0.00723 68 44.491  <.0001 
##  High - Low       0.455 0.00653 68 69.754  <.0001 
##  Medium - Low     0.134 0.00728 68 18.355  <.0001 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

3.6 Growth rate and He

#Best model with genetic diversity
mod0 <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population), 
                   data=data_phenotyping)

mod0 <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_remain_end),])
mod1 <- lme4::lmer(Lambda ~ He_remain_end + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_remain_end),])
mod2 <- lme4::lmer(Lambda ~ Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_remain_end),])
mod3 <- lme4::lmer(Lambda ~ log(He_remain_end) + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_remain_end),])
mod4 <- lme4::lmer(Lambda ~ 1 + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_remain_end),])
mod5 <- lme4::lmer(Lambda ~ exp(He_remain_end) + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_remain_end),])
MuMIn::model.sel(mod0, mod1,mod2,mod3,mod4,mod5)
## Model selection table 
##        (Int) Blc Gnt_Dvr He_rmn_end log(He_rmn_end) exp(He_rmn_end)
## mod3  2.6620   +                               1.77                
## mod1  0.2305   +              2.471                                
## mod5 -0.4063   +                                              1.151
## mod0  2.6630   +       +                                           
## mod2  1.9410   +                                                   
## mod4  1.9410   +                                                   
##                  family df   logLik  AICc delta weight
## mod3 gaussian(identity)  6 -283.053 578.5  0.00  0.473
## mod1 gaussian(identity)  6 -283.276 579.0  0.45  0.379
## mod5 gaussian(identity)  6 -284.363 581.2  2.62  0.128
## mod0 gaussian(identity)  7 -285.136 584.8  6.31  0.020
## mod2 gaussian(identity)  5 -297.218 604.7 26.21  0.000
## mod4 gaussian(identity)  5 -297.218 604.7 26.21  0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
## '1 | Population'
#Makes sense: best model log He
x <- seq(1:100)
y <- seq(1001:1100)
plot(log(x),y)

PLOT_He_withextinct

#############################################################
##################        Analysis        ###################
#############################################################

mod3 <- lme4::lmer(Lambda ~ log(He_remain_end) + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_remain_end),])
mod4 <- lme4::lmer(Lambda ~ 1 + Block + (1|Population), 
                   data=data_phenotyping[!is.na(data_phenotyping$He_remain_end),])
anova(mod3,mod4,test="Chisq") # difference
## Data: data_phenotyping[!is.na(data_phenotyping$He_remain_end), ]
## Models:
## mod4: Lambda ~ 1 + Block + (1 | Population)
## mod3: Lambda ~ log(He_remain_end) + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## mod4    5 599.98 616.53 -294.99   589.98                         
## mod3    6 571.83 591.68 -279.91   559.83 30.158  1  3.983e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Lambda ~ log(He_remain_end) + Block + (1 | Population)
##    Data: data_phenotyping[!is.na(data_phenotyping$He_remain_end), ]
## 
## REML criterion at convergence: 566.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4098 -0.6488 -0.0811  0.5245  3.3605 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Population (Intercept) 0.2616   0.5115  
##  Residual               0.7609   0.8723  
## Number of obs: 202, groups:  Population, 73
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)          2.6623     0.1900  14.015
## log(He_remain_end)   1.7700     0.3028   5.846
## BlockBlock4          0.2491     0.2060   1.209
## BlockBlock5         -0.1118     0.2391  -0.468
## 
## Correlation of Fixed Effects:
##             (Intr) l(H__) BlckB4
## lg(H_rmn_n)  0.634              
## BlockBlock4 -0.617 -0.103       
## BlockBlock5 -0.569 -0.148  0.454
#############################################################
##################        Predict         ###################
#############################################################

data_predict_lambda <- data.frame(He_remain_end =
  seq(min(data_phenotyping[!is.na(data_phenotyping$He_remain_end),]$He_remain_end),
      max(data_phenotyping[!is.na(data_phenotyping$He_remain_end),]$He_remain_end), 
      0.01),
                                   Block = "Block4")

data_predict_lambda$lambda_predict <- predict(mod3, 
                                        type = "response",
                                        re.form = NA,
                                        newdata = data_predict_lambda)
 
   

PLOT_He_expect <- ggplot2::ggplot() +
  geom_line(data = data_predict_lambda, 
                                aes(x = He_remain_end, y = lambda_predict), 
                                linetype = "longdash", col = "grey40", size = 1.25) +
  geom_point(data = SUM_POP_He_ALL_extinct, 
                               aes(x = He_remain_end, y = Lambda, 
                         colour = Genetic_Diversity, 
                         shape = Extinction_finalstatus, 
                         size = N), alpha = 0.7) +
  ylab("Growth rate") +
  xlab("Remaining heterozygosity") +
  scale_shape_manual(values = c(19, 1)) + 
  theme_LO + 
  labs(color="Genetic diversity") +
  labs(shape="Extinct population") +
  labs(size="Replicates") 

PLOT_He_expect

# 
# cowplot::save_plot(here::here("figures","Fitness_He_final.pdf"),   PLOT_He_withextinct,
#           base_height = 12/cm(1), base_width = 18/cm(1), dpi = 610)

3.7 Adults G6

m0 <- lme4::glmer(Adults ~ Genetic_Diversity  + (1|Population) + (1|ID_Rep:Population),
                family = "poisson", data = data_phenotyping)
summary(m0)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
##    Data: data_phenotyping
## 
##      AIC      BIC   logLik deviance df.resid 
##   2154.4   2171.3  -1072.2   2144.4      212 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1441 -0.1230  0.0254  0.1360  0.3923 
## 
## Random effects:
##  Groups            Name        Variance Std.Dev.
##  ID_Rep:Population (Intercept) 0.2041   0.4518  
##  Population        (Intercept) 0.2807   0.5298  
## Number of obs: 217, groups:  ID_Rep:Population, 217; Population, 77
## 
## Fixed effects:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               4.2792     0.1105  38.715  < 2e-16 ***
## Genetic_DiversityMedium  -0.7386     0.1812  -4.076 4.59e-05 ***
## Genetic_DiversityLow     -0.5250     0.1664  -3.156   0.0016 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Gnt_DM
## Gntc_DvrstM -0.611       
## Gntc_DvrstL -0.667  0.413
# Equivalent to:
# m0 <- lme4::glmer(Adults ~ Genetic_Diversity  + (1|Population/ID_Rep),
#                 family = "poisson", data = data_5week)


#dispersion_lme4::glmer(m0) #dispersion ok 
#Note: (1|School) + (1|Class:School) is equivalent to (1|School/Class)

m0 <- lme4::glmer(Adults ~ Genetic_Diversity  + (1|Population) + (1|ID_Rep:Population),
                family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1  + (1|Population) + (1|ID_Rep:Population),
                family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") # difference
## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Population) + (1 | ID_Rep:Population)
## m0: Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
##    npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)    
## m1    3 2167.5 2177.6 -1080.7   2161.5                        
## m0    5 2154.4 2171.3 -1072.2   2144.4 17.04  2  0.0001994 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m0)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
##    Data: data_phenotyping
## 
##      AIC      BIC   logLik deviance df.resid 
##   2154.4   2171.3  -1072.2   2144.4      212 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1441 -0.1230  0.0254  0.1360  0.3923 
## 
## Random effects:
##  Groups            Name        Variance Std.Dev.
##  ID_Rep:Population (Intercept) 0.2041   0.4518  
##  Population        (Intercept) 0.2807   0.5298  
## Number of obs: 217, groups:  ID_Rep:Population, 217; Population, 77
## 
## Fixed effects:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               4.2792     0.1105  38.715  < 2e-16 ***
## Genetic_DiversityMedium  -0.7386     0.1812  -4.076 4.59e-05 ***
## Genetic_DiversityLow     -0.5250     0.1664  -3.156   0.0016 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Gnt_DM
## Gntc_DvrstM -0.611       
## Gntc_DvrstL -0.667  0.413
############ 
###### Posthoc
############ 
emmeans::emmeans(m0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean    SE  df asymp.LCL asymp.UCL
##  High                4.28 0.111 Inf      4.02      4.54
##  Medium              3.54 0.143 Inf      3.20      3.88
##  Low                 3.75 0.124 Inf      3.46      4.05
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE  df z.ratio p.value
##  High - Medium    0.739 0.181 Inf  4.076  0.0001 
##  High - Low       0.525 0.166 Inf  3.156  0.0046 
##  Medium - Low    -0.214 0.189 Inf -1.132  0.4943 
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates
#test double nested
m0 <- lme4::glmer(Adults ~ Genetic_Diversity  + (1|Genetic_Diversity/Population/ID_Rep),
                family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1  + (1|Genetic_Diversity/Population/ID_Rep),
                family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") 
## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Genetic_Diversity/Population/ID_Rep)
## m0: Adults ~ Genetic_Diversity + (1 | Genetic_Diversity/Population/ID_Rep)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## m1    4 2160.9 2174.4 -1076.4   2152.9                       
## m0    6 2156.4 2176.7 -1072.2   2144.4 8.4441  2    0.01467 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m0 <- lme4::glmer(Adults ~ Genetic_Diversity  + (1|Genetic_Diversity/ID_Rep),
                family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1  + (1|Genetic_Diversity/ID_Rep),
                family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") 
## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Genetic_Diversity/ID_Rep)
## m0: Adults ~ Genetic_Diversity + (1 | Genetic_Diversity/ID_Rep)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## m1    3 2189.9 2200.0 -1091.9   2183.9                        
## m0    5 2183.0 2199.9 -1086.5   2173.0 10.858  2   0.004387 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.8 Adults G2

m0 <- glm(Nb_adults ~ Genetic_Diversity + Block,  family = "poisson", data = data_G2)
summary(m0)
## 
## Call:
## glm(formula = Nb_adults ~ Genetic_Diversity + Block, family = "poisson", 
##     data = data_G2)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -18.7057   -2.7621    0.2126    3.3983   11.4454  
## 
## Coefficients:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              5.73890    0.01497 383.312  < 2e-16 ***
## Genetic_DiversityMedium -0.19408    0.01628 -11.920  < 2e-16 ***
## Genetic_DiversityLow    -0.19278    0.01460 -13.205  < 2e-16 ***
## BlockBlock4              0.10767    0.01579   6.817  9.3e-12 ***
## BlockBlock5              0.17743    0.01600  11.088  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2379.6  on 83  degrees of freedom
## Residual deviance: 2027.9  on 79  degrees of freedom
## AIC: 2667.2
## 
## Number of Fisher Scoring iterations: 4
m1<- glm(Nb_adults ~ Block,  family = "poisson", data = data_G2)
anova(m0,m1,test="Chisq") # difference
## Analysis of Deviance Table
## 
## Model 1: Nb_adults ~ Genetic_Diversity + Block
## Model 2: Nb_adults ~ Block
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1        79     2027.9                          
## 2        81     2241.4 -2  -213.52 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############ 
###### Posthoc
############ 
emmeans::emmeans(m0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean      SE  df asymp.LCL asymp.UCL
##  High               5.834 0.01051 Inf     5.809     5.859
##  Medium             5.640 0.01243 Inf     5.610     5.670
##  Low                5.641 0.01022 Inf     5.617     5.666
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate     SE  df z.ratio p.value
##  High - Medium   0.1941 0.0163 Inf 11.920  <.0001 
##  High - Low      0.1928 0.0146 Inf 13.205  <.0001 
##  Medium - Low   -0.0013 0.0161 Inf -0.081  0.9964 
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates

3.9 Evolution of growthrate over time

############ 
###### Distribution
############
hist(data_evolution_Lambda$Lambda, breaks=50)

# Square root
msqrt <- lme4::lmer(sqrt(Lambda) ~ Genetic_Diversity*Phenotyping  + 
                Block + (1|Population), 
              data = data_evolution_Lambda)
summary(msqrt)
## Linear mixed model fit by REML ['lmerMod']
## Formula: sqrt(Lambda) ~ Genetic_Diversity * Phenotyping + Block + (1 |  
##     Population)
##    Data: data_evolution_Lambda
## 
## REML criterion at convergence: 252.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4979 -0.4737  0.0672  0.5620  2.4161 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Population (Intercept) 0.0295   0.1718  
##  Residual               0.1040   0.3225  
## Number of obs: 301, groups:  Population, 88
## 
## Fixed effects:
##                                           Estimate Std. Error t value
## (Intercept)                               1.779696   0.081683  21.788
## Genetic_DiversityMedium                  -0.171161   0.103522  -1.653
## Genetic_DiversityLow                     -0.150003   0.094035  -1.595
## PhenotypingFinal                         -0.243939   0.070714  -3.450
## BlockBlock4                               0.102764   0.063883   1.609
## BlockBlock5                               0.003821   0.070206   0.054
## Genetic_DiversityMedium:PhenotypingFinal -0.251050   0.110515  -2.272
## Genetic_DiversityLow:PhenotypingFinal    -0.154999   0.101256  -1.531
## 
## Correlation of Fixed Effects:
##             (Intr) Gnt_DM Gnt_DL PhntyF BlckB4 BlckB5 G_DM:P
## Gntc_DvrstM -0.600                                          
## Gntc_DvrstL -0.655  0.506                                   
## PhntypngFnl -0.685  0.534  0.588                            
## BlockBlock4 -0.462  0.054  0.041  0.028                     
## BlockBlock5 -0.424  0.006  0.009  0.002  0.488              
## Gntc_DvM:PF  0.428 -0.745 -0.377 -0.641 -0.028  0.050       
## Gntc_DvL:PF  0.444 -0.373 -0.733 -0.699  0.012  0.088  0.453
# log
mlog <- lme4::lmer(log(Lambda+1) ~ Genetic_Diversity*Phenotyping  +
               Block + (1|Population),
              data=data_evolution_Lambda)

# Normal
mnorm <- lme4::lmer(Lambda ~ Genetic_Diversity*Phenotyping  + 
                Block + (1|Population),
              data=data_evolution_Lambda)

## Compare AIC
AIC(mlog, msqrt,mnorm)
##       df      AIC
## mlog  10 207.0029
## msqrt 10 272.0671
## mnorm 10 848.9685
#Square root a better fits to the data 
#But we compare different values: Growth rate and Growth rate+1

#Distribution Residuals
hist(residuals(mlog),col="yellow",freq=F)

hist(residuals(msqrt),col="yellow",freq=F)  #Seems the better fit

hist(residuals(mnorm),col="yellow",freq=F)

#Test normality
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(mlog), "norm"))
g1$kstest 
## 1-mle-norm 
## "rejected"
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(msqrt), "norm"))
g1$kstest 
##     1-mle-norm 
## "not rejected"
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(mnorm), "norm"))
g1$kstest
##     1-mle-norm 
## "not rejected"
# QQ plot
# plot(mlog, which = 2) 
# plot(msqrt, which = 2) #Seems the better fit!
# plot(mnorm, which = 2) #Seems the better fit but very high AIC 
qqnorm(resid(mlog))
qqline(resid(mlog))

qqnorm(resid(msqrt))
qqline(resid(msqrt))

qqnorm(resid(mnorm))
qqline(resid(mnorm)) #best fit

############ 
###### Model selection
############

m0 <- lme4::lmer(Lambda ~ Genetic_Diversity*Phenotyping  + 
                Block + (1|Population) ,
              data=data_evolution_Lambda)


m1 <- lme4::lmer(Lambda ~ Genetic_Diversity+Phenotyping  + 
                Block + (1|Population) ,
              data=data_evolution_Lambda)

anova(m0,m1,test="Chisq") #Remove interaction
## Data: data_evolution_Lambda
## Models:
## m1: Lambda ~ Genetic_Diversity + Phenotyping + Block + (1 | Population)
## m0: Lambda ~ Genetic_Diversity * Phenotyping + Block + (1 | Population)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m1    8 832.93 862.59 -408.46   816.93                     
## m0   10 834.25 871.32 -407.12   814.25 2.6798  2     0.2619
#Same result with: 
# m1 <- lmer(Lambda ~ Genetic_Diversity+Phenotyping  + 
#                 (1|Block) + (1|Population) + (1|ID_Rep:Population),
#               data=data_evolution_Lambda)
# OR
# m1 <- lmer(Lambda ~ Genetic_Diversity+Phenotyping  + 
#                (1|Block) + (1|Population) ,
#               data=data_evolution_Lambda)

m2 <- lme4::lmer(Lambda ~ Phenotyping  + 
                Block + (1|Population),
              data=data_evolution_Lambda)
m3 <- lme4::lmer(Lambda ~ Genetic_Diversity  + 
                Block + (1|Population),
              data=data_evolution_Lambda)

anova(m1,m2,test="Chisq") #Keep phenotyping
## Data: data_evolution_Lambda
## Models:
## m2: Lambda ~ Phenotyping + Block + (1 | Population)
## m1: Lambda ~ Genetic_Diversity + Phenotyping + Block + (1 | Population)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## m2    6 855.30 877.54 -421.65   843.30                         
## m1    8 832.93 862.59 -408.46   816.93 26.369  2   1.88e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m1,m3,test="Chisq") #Keep genetic diversity
## Data: data_evolution_Lambda
## Models:
## m3: Lambda ~ Genetic_Diversity + Block + (1 | Population)
## m1: Lambda ~ Genetic_Diversity + Phenotyping + Block + (1 | Population)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## m3    7 896.49 922.44 -441.24   882.49                         
## m1    8 832.93 862.59 -408.46   816.93 65.558  1  5.642e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############ 
###### Posthoc
############ 
mod_final <- lme4::lmer(Lambda ~ Genetic_Diversity+Phenotyping  + 
                Block + (1|Population),
              data=data_evolution_Lambda)

emmeans::emmeans(mod_final, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean    SE   df lower.CL upper.CL
##  High                3.05 0.120 69.4     2.76     3.34
##  Medium              2.14 0.143 87.6     1.79     2.48
##  Low                 2.32 0.123 99.8     2.02     2.62
## 
## Results are averaged over the levels of: Phenotyping, Block 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE   df t.ratio p.value
##  High - Medium    0.916 0.184 72.8  4.988  <.0001 
##  High - Low       0.734 0.170 78.2  4.316  0.0001 
##  Medium - Low    -0.182 0.187 87.5 -0.972  0.5963 
## 
## Results are averaged over the levels of: Phenotyping, Block 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates
emmeans::emmeans(mod_final, list(pairwise ~ Phenotyping), adjust = "tukey") #
## $`emmeans of Phenotyping`
##  Phenotyping emmean     SE  df lower.CL upper.CL
##  Initial       3.00 0.1072 254     2.76     3.24
##  Final         2.01 0.0837 109     1.82     2.20
## 
## Results are averaged over the levels of: Genetic_Diversity, Block 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 2 estimates 
## 
## $`pairwise differences of Phenotyping`
##  contrast        estimate    SE  df t.ratio p.value
##  Initial - Final    0.993 0.117 250 8.474   <.0001 
## 
## Results are averaged over the levels of: Genetic_Diversity, Block 
## Degrees-of-freedom method: kenward-roger
#